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^\ ' Interferometry at radio frequencies between Earth-based receivers separated by intercontinental 

distances has made significant contributions to astrometry and geophysics during the past three 
decades. Analyses of such Very Long Baseline Interferometric (VLBI) experiments now permit 
measurements of relative positions of points on the Earth's surface, and angles between celestial 

Q objects, at the levels of better than 1 cm and 1 nanoradian, respectively. The relative angu- 
lar positions of extragalactic radio sources inferred from this technique presently form the best 
realization of an inertial reference frame. This review summarizes the current status of radio 
interferometric measurements for astrometric and geodetic applications. It emphasizes the theo- 
retical models that are required to extract results from the VLBI observables at present accuracy 
levels. An unusually broad cross-section of physics contributes to the required modeling. Both 
special and general relativity need to be considered in properly formulating the geometric part 
^ ' of the propagation delay. While high-altitude atmospheric charged particle (ionospheric) effects 

oo : are easily calibrated for measurements employing two well separated frequencies, the contribution 

CO ■ of the neutral atmosphere at lower altitudes is more difficult to remove. In fact, mismodeling 

' of the troposphere remains the dominant error source. Plate tectonic motions of the observing 

^Nj ' stations need to be taken into account, as well as the non-point-like intensity distributions of many 

sources. Numerous small periodic and quasi-periodic tidal effects also make important contribu- 
^ tions to space geodetic observables at the centimeter level, and some of these are just beginning 

I to be characterized. Another area of current rapid advances is the specification of the orientation 

1 of the Earth's spin axis in inertial space: nutation and precession. Highlights of the achievements 

' of VLBI are presented in four areas: reference frames. Earth orientation, atmospheric effects on 

^1^' microwave propagation, and relativity. The order-of-magnitude improvement of accuracy that was 

achieved during the last decade has provided essential input to geophysical models of the Earth's 
internal structure. Most aspects of VLBI modeling are also directly applicable to interpretation of 
other space geodetic measurements, such as active and passive ranging to Earth-orbiting satellites, 
' interplanetary spacecraft, and the Moon. 
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Symbols and abbreviations 

Aij nutation amplitudes in longitude 

B baseline vector 

Bij nutation amplitudes in obliquity 

BIH Bureau International de I'Heure 

BIPM Bureau International des Poids et Mesures 

CDF Crustal Dynamics Project 

CIO Conventional International Origin 

Cjpj planetary nutation amplitudes in longitude (Eq. 3.116) 

Cgj planetary nutation amplitudes in obliquity (Eq. 3.117) 

c speed of light 

D mean elongation of Moon from Sun 

E elevation angle 

F latitude argument of Moon 

/ flattening factor of the Earth 

/ pressure loading factor 

G universal gravitational constant 

GPS Global Positioning System 

GSFC Goddard Space Flight Center 

g gravitational acceleration at Earth's surface 

g angle for barycentric dynamic time (Eq. (3.32) 

H hour angle 

h altitude above reference ellipsoid 

hr hour angle of mean equinox of date 
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hi vertical Love number [i = 2,3: quadrupole, octupole) 

hi^2 ionosphere or troposphere height Umits 

lAU International Astronomical Union 

ICRF International Celestial Reference Frame 

lERS International Earth Rotation Service 

IRIS International Radio Interferometric Surveying 

ITRF International Terrestrial Reference Frame 

lUGG International Union of Geodesy and Geophysics 

JD Julian date 

k unit vector in signal propagation direction 

LLR Lunar Laser Ranging 

I mean anomaly of Moon 

I' mean anomaly of Sun 

/e mean longitude of Earth 

I J mean longitude of Jupiter 

Im mean longitude of Mars 

Zs mean longitude of Saturn 

ly mean longitude of Venus 

k horizontal Love number {i = 2,3: quadrupole, octupole) 

mas milliarcsecond 

MERIT Monitor Earth Rotation and Intercompare Techniques 

MIT Massachusetts Institute of Tcichnology 

m speed of precession in right ascension 

nip mass of body p 

M„ first moment of wet troposphere refractivity 

A^d,w Mapping functions for dry, wet troposphere 

N nutation transformation matrix 

N atmospheric refractivity 

N' lunar node argument 

NASA National Aeronautics and Space Administration 

NEOS National Earth Orientation Service 

NMF Niell mapping function 

NNR No Net Rotation 

NOAA National Oceanic and Atmospheric Administration 

NRAO National Radio Astronomy Observatory 

Nuvel New velocity model for plate tectonics 

n refractive index 

n speed of precession in declination 

P precession transformation matrix 

p extended pressure anomaly 

Pa general precession 

Phs lunisolar precession 

PPL planetary precession 

PSTD reference atmospheric pressure 

Q rotation matrix for terrestrial to celestial transformation 

RA Right ascension 

Re Earth equatorial radius 

Rc Earth center coordinates in SSB frame 

Rp position of perturbing body in geocentric celestial system 

Rp geocentric unit vector to perturbing body 

i?EG distance from Earth to gravitating body G 

r unit vector in radial direction 

ro classical electron radius 

To station position in terrestrial system, excluding motions 

Tc station position in geocentric celestial system 
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phase shifted station position 

geocentric unit vector to station position 

station position in terrestrial frame 

position of station 1,2 in SSB frame 

station position of date 

station cylindrical radius from spin axis 

slant range factor 

S-band + X-band 

Space Geodesy Project 

Satellite Laser Ranging 

Solar System Barycenter 

planetary nutation amplitudes in longitude 

planetary nutation amplitudes in obliquity 

geocentric celestial unit vector to source 

SSB celestial unit vector to source 

Temps Dynamique Barycentrique 

Terrestrial Dynamic Time 

time in centuries since J2000.0 

reference time 

time of arrival of wave front at station 1 

proper time of arrival of wave front at station 1 

time of arrival of wave front at station 2 

time of emission by source 

light transit time 

UTl transformation matrix 

gravitational potential 

Universal Time 1 

Universal Time Coordinated 

Universal Time and Polar Motion 

projections of B on plane of sky 

transformation matrix from local to equatorial frame (Eqs. 3.61) 

Vertical, East, North (local geodetic coordinates) 

Very Large Array 

Very Long Baseline Array 

Very Long Baseline Interferometry 

astronomical argument of tidal constituent i 

atmospheric temperature lapse rate 

Water Vapor Radiometer 

polar motion transformation matrix, x component 

Cartesian coordinates of station i 

Cartesian coordinates of station i at reference time to 

Cartesian velocities of station i 

polar motion transformation matrix, y component 

auxiliary angle for precession (Eq. 3.129) 

zenith dry, wet troposphere delay 

Zhu, Mathews, Oceans, Anelasticity (nutation model) 

right ascension 

equation of the equinoxes 

velocity of station 1, 2 

general relativity (Parametrized Post-Newtonian) gamma factor 

total tidal shift in terrestrial coordinate system 

atmospheric loading station position shift 

ionosphere contribution to group delay 

height of antenna reference point above surface 

ocean loading station position shift 
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ij.) phase shift of sohd tidal effect 

fl mean longitude of ascending lunar node 

n perturbation transformation matrix 

lue rotational speed of Earth 

LV[ free core nutation frequency 

col Cartesian components of angular velocity of tectonic plate j 



I. INTRODUCTION 

Astrometry and geodesy have undergone a revolution during the past three decades. This revolution was initiated 
in the mid-1960s by the advent of interferometry at radio frequencies using antennas separated by thousands of 
kilometers. Developments during the following decades refined this "Very Long Baseline Interferometry" (VLBI) 
technique to reach its present capability of point positioning on the Earth's surface at the centimeter level, and 
angular positioning for point-like extragalactic radio sources at the sub-milliarcsecond (nanoradian) level. These 
refinements are expected to continue until insurmountable problems are reached, probably in the sub-millimeter 
accuracy regime. Geodetic measurements on the Earth's surface have been enormously expanded and densified by 
satellite tracking methods developed during the past decade, but VLBI remains the prime technique for astrometry 
employing natural radio sources. It is unique in its ability to measure the Earth's orientation in an inertial frame of 
reference. 

Detection and processing of the extremely weak signals from distant radio sources requires long integration times, 
sophisticated data acquisition systems, and customized computers. Analyses of VLBI observables involve an uncom- 
monly broad cross-section of the disciplines of physical science, ranging from consideration of the effects of the Earth's 
internal structure on its dynamics, to the tectonic plate motions and various terrestrial tidal effects, to quantification 
of turbulence in the atmosphere, to special relativistic description of the radio signals traveling from the distant 
sources, to general relativistic retardation. With the exception of gravimetric and oceanographic experiments, VLBI 
is perhaps the most demanding technique for many aspects of global Earth models. It therefore provides results that 
are used in a number of related fields. While extension of experiments to platforms in Earth orbit and elsewhere in 
the Solar System is in its infancy (Burke, 1991; Hirabayashi et aL, 1991; Adam, 1995), the present review is limited 
to Earth-based observations. A large portion of the modeling we shall discuss here is common to VLBI and Earth- 
satellite techniques such as laser ranging to reflectors on the Moon (lunar laser ranging, LLR) and ranging to passive 
(Lagcos, SLR) and active satellites (Global Positioning System, GPS). For these reasons, it is important to review 
the complete details of current models of VLBI observables. Following the principle that the accuracy of a theoretical 
model should exceed the accuracy of the experiments which it interprets by at least an order of magnitude, the VLBI 
model should ideally be complete at the 1 pic;os{X"ond lcvc;l. Unfortunatc;ly this requirenic;nt is currently not satisfied 
by a number of parts of the model. As experimental techniques are refined, possibly to improve the accuracy by 
another order of magnitude, numerous aspects of the model will need to be re-examined. 

The field of radio astronomy was initiated in the 1930s by the discovery of celestial radio emission by Jansky 
(1932; 1933). Development of instrumentation continued, and in the 1940s Reber (1940; 1944) used a parabolic 
antenna to construct the first maps of radio emission from the Milky Way galaxy at frequencies of 160 and 480 MHz. 
Subsequent evolution of radio astronomy through the post-war years has been traced in the comprehensive collection 
of memoirs edited by Sullivan (1984). Following the discovery of the powerful non-thermal radio sources called 
quasars (Schmidt, 1963), which are extremely distant and emit vast quantities of energy, intense interest developed 
in characterizing their intrinsic nature, as well as in using them as the basis for a novel astronomical reference frame. 
The 3rd Cambridge survey (Edge et al., 1959) and the Parkes survey (Bolton et al, 1979) identified many of the radio 
sources used today for geodetic and astrometric VLBI. 

The basic ideas of the VLBI technique were first demonstrated in experiments detecting decamctric bursts of Jupiter 
by a group at the University of Florida in 1965 (Carr et al, 1965; May and Carr, 1967). Subsequent history can be 
traced from the first extragalactic interferometric observations by the NRAO-Arecibo group (Bare et al., 1967) at 
0.6 GHz on a 220-km baseline, the MIT-NRAO group (Moran et al. (1967); 1.7 GHz, 845 km), and the Canadian 
Long Baseline Interferometry group (Broten et al. (1967); 0.4 GHz, 3074 km). Early work in geodesy, astrometry, 
and clock synchronization was done in 1969 (Hinteregger et al., 1972; Cohen, 1972), yielding accuracies in distances 
of 2—5 meters and in source positions on the order of 1 arc second. After the first decade spent in developing 
the technique, numerous applications matured in succeeding years using networks of fixed antennas. Extension to 
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transportable antennas was also implemented during the 1970s (MacDoran, 1974; Ong et ai, 1976; Davidson and 
Trask, 1985). During that period, several basic reviews of VLBI instrumentation, experiments, and analysis appeared 
as chapters in the series "Methods of Experimental Physics: Astrophysics" (Meeks, 1976). Some Ph. D. dissertations 
of that era arc those of Whitney (1974), Robertson (1975), and Ma (1978). After the currently standard Mark III 
VLBI data acquisition system was developed (Rogers et ai, 1983), large-scale international cooperation began in 1979 
with the estabhshment of the Crustal Dynamics Project (CDP) by NASA (Bosworth et al, 1993). In the 1980s, 
radio interferometry progressed to the point of yielding some of the most accurately measured physical quantities. 
Increasingly larger and more sensitive antennas were used to extend observations to ever weaker sources of emission, 
and the results have made important contributions to our understanding of a wide variety of astrophysical and 
geodynamic phenomena. Such advances in instrumentation culminated in the construction of instruments such as the 
Very Large Array (VLA) in 1980 (Hjellming and Bignell, 1982), and the Very Long Baseline Array (VLB A) in 1993 
(Napier, 1994). 

Useful references which provide an introduction to the principles of very long baseline interferometry are 
several articles in the above-mentioned volume 12 of "Methods of Experimental Physics" (Rogers, 1976; Poo- 
ley, 1976; Moran, 1976; Vessot, 1976; Shapiro, 1976), the book by Thompson, Moran and Swenson (1986), and 
two reports by Thomas (1981; 1987). Additional background material can be foimd in the review of geophysical 
VLBI applications by Robertson (1991), and in more detail and more recently, in the three-volume compilation that 
summarizes the accomplishments of NASA's Crustal Dynamics Project during the past decade (Smith and Tur- 
cotte, 1993). Various periodic reviews, including the quadrennial lUGG geodesy reports (Clark (1979), Carter (1983); 
Robertson (1987); Ray (1991); Herring (1995)), the lERS annual reports {e.g. lERS, 1996a), and the yearly reports 
on CDP results [e.g. Ma et al. (1992), more recently at http://lupus.gsfc.nasa.gov/vlbi.litml], are also good 
sources for detailed information concerning current VLBI techniques and results. The International Earth Rotation 
Service (lERS) periodically publishes a compilation of standard models recommended for analyses of space geodetic 
data (1992; 1996b); also available at http://inaia.usno.navy.niil/conventions.htinl. In large part, the model 
description in the present paper is in agreement with the specifications of these "lERS Standards" . 

The intent of this paper is to review astrometric and geodetic VLBI from a description of the experiments through 
analysis of the results, with emphasis on the details of theoretical modeling. Section II summarizes experiment design, 
single-station signal collection and recording, signal combination, and observable extraction. This section also reviews 
the general approach to extracting parameters of astrometric and geodetic interest from the observables. The three 
succeeding sections (III-V) consider the theoretical model in detail, partitioning the interferometric delay model into 
three major components: geometry, instrumentation (clocks), and atmosphere, and presenting the best current models 
of each of these components. The longest section (III) deals with the purely geometric portion of the delay, wave 
front curvature, and gravitational bending, and considers time definitions, tectonic motion, tidal and source structure 
effects, coordinate frames. Earth orientation (universal time and polar motion), nutation, precession. Earth orbital 
motion, and antenna offsets. Sections IV and V discuss the non-gcomctric components of the model: clocks and 
the atmosphere. Section VI presents a brief overview of the accomplishments of VLBI and its contributions to our 
understanding of astrometry and geophysics during the past two decades. Section VII outlines model improvements 
that we anticipate will be required by more accurate data in the future. 

The most convenient units for use in astrometric and geodetic VLBI are millimeters, picoseconds, and nanoradians. 
The first two are connected by the speed of light = 299,792,458 m/s «1 mm/3 ps, and will be used interchangeably. 
Most commonly for the purposes of illustration, a 10,000-km baseline will be considered. A length change of 1 cm on 
such a baseline is equivalent to an angular change of 1 nrad (1 part per billion, ppb) or 0.2 milliarcseconds (mas); the 
equivalent conversion factors in millimeters are 10 mm/nrad or 50 mm/mas. 

II. EXPERIMENTS AND ANALYSIS 

Interferometric measurements permit resolution on the order of the wavelength of the radiation. For radio waves 
at the typical astrometric/geodetic frequencies of 2.3 GHz (S band) and 8.4 GHz (X band), the wavelengths are 
13 and 3.5 cm, respectively. Some recent activity is focusing on even higher frequencies {e.g. Lebach et al., 1995), 
and numerous interferometric observations are planned at millimeter wavelengths, spanning the gap between the 
radio and optical regions. For experiments that are dominated by random errors, gains can be achieved by repeated 
measurements, thereby reducing the errors substantially. Thus, the accuracy at X band may be expected to surpass 
1 cm in many cases. 

Early radio interferometric experiments employed antennas that were separated by distances on the order of 1 
km. These experiments had the advantage that the stations shared a frequency standard, and signal correlation 
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was performed in real time {e.g. Batty et al, 1982). Such an arrangement is now known as a "connected element 
interferometer" to distinguish it from interferometers with separated elements that are not in direct communication. 
A connected element interferometer is a close analogue of the Michelson stellar interferometer which manipulates 
signals with mirrors to produce a physical interference pattern at the detector. Frequency standards ("clocks") 
and the technology associated with measurements at centimeter wavelengths improved dramatically as a result of 
radar development during World War II. It became possible to separate the two components of the interferometer 
with independent clocks and instrumentation, and to use high-speed data recording techniques, thus eliminating the 
need for close physical proximity. With such an arrangement there are no physical real-time interference "fringes" , 
although the terminology of optical interferometry has been carried over to the radio wavelength region. Remote radio 
interferometry is now commonplace. With the two antennas separated by long baselines, the experimental observables 
are no longer obtained in real time, but are instead generated by subsequent analysis of recorded information. Such 
experiments have become known as "very long" baseline interferometry or VLBI. With the recent explosive advances 
in communications, computing, and data storage technology, connected element experiments over intercontinental 
distances may again become standard in the future. 

Large dish antennas (20 to 100 m in diameter) that are used for radio astronomy and VLBI are found in approxi- 
mately 50 globally distributed locations. An equivalent number of sites have been occupied by smaller transportable 
antennas during the past two decades. The premier sites are at relatively high altitudes and in arid climates, which 
were selected to minimize the impact of atmospheric turbulence on measurement quality. 



A. Experiment design and data collection 



A diagram of a contemporary VLBI experiment is sketched in Fig. 1. Two antennas, separated by a baseline B, are 
simultaneously pointed at the same radio source, and detect the incoming wave front propagating along unit vector k. 





rrrrr 



CORRELATOR 



FIG. 1. Schematic diagram of a VLBI experiment. Waveforms are shown impinging from the direction k on two antennas 

separated by a baseline B on the Earth's surface. They are followed through the data acquisition system to the point where 
the correlator determines the delay t. The signal waveforms are exaggerated for effect. The actual waveforms are random 
Gaussian processes. 
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More than two antennas (stations) commonly participate in an experiment, but its fundamentals can be adequately 
illustrated by considering only one of the many pairs of stations. Prior to its arrival at the first antenna, the signal 
has already been affected by the electromagnetic and gravitational fields that it encountered along its journey through 
interstellar space, the Solar System, and the Earth's atmosphere. Between this time and the time of reception at 
the second antenna, any motions of the Earth also make their influence felt. All these varied aspects need to be 
taken into account in analyses of the experiment. We note that the VLBI technique is differential in the sense that 
the observables depend on the relaMve station locations and source positions, but they are fairly insensitive to the 
stations' absolute geocentric locations and the origin of right ascension (celestial longitude). Other techniques {e.g., 
SLR, GPS) must supply such information, which minimally consists of the geocentric position of one of the observing 
stations at a given epoch, and the right ascension of one source. 

The information content of the VLBI observables depends critically on the design of the experiment {e.g. Shapiro, 
1976). The observing schedule (station network, observation epochs, and local azimuth and elevation angles of the 
selected sources) plays a crucial role in determining the types and precision of parameters that can be extracted in 
the subsequent data analyses. At a bare minimum there must be at least one observable for each parameter to be 
estimated. The experimenter must first decide how many parameters he wishes to measure. To consider a simple 
example: after having fixed the station location and clock parameters of a reference station, a typical astrometric or 
geodetic analysis will want to estimate, at a minimum, for each of the remaining stations three site coordinates, a 
zenith tropospheric delay parameter, and linear clock parameters (offset and rate). In addition to these station-specific 
parameters, a number of "global" parameters are also usually estimated. These are common to the entire network, and 
describe the orientation of the Earth (5 parameters) and the positions of the sources (2 angular parameters per source); 
the right ascension of one source must be fixed to establish the origin of the celestial reference frame. Thus for an 
experiment with N stations and M sources, there would be Np parameters to determine, where Np = 6N + 2M — 1. 
In recent years, analyses have tended to estimate increasingly larger numbers of station-specific parameters for the 
clocks and tropospheres, thereby increasing the coefficient of the TV term many-fold. 

Note that the observable delay in arrival times is approximately proportional to the scalar product of the baseline and 
signal propagation vectors, and thus contains information concerning both the station and source coordinates. With 
sufficiently redundant measurements, therefore, both the station and source coordinates are accessible to estimation, 
subject to the limitations imposed by only observing the scalar product. Naturally, it is highly desirable to have more 
observables than the minimum {Np) required for a non-singular solution in the analysis. The additional observables 
permit reduction of the formal statistical uncertainties in the estimated parameters. They also allow the level of 
systematic errors to be characterized using various tests of the repeatability of the results. Such tests may include 
root-mean-square of the residuals (observed minus theoretical) and repeatability of estimated parameters derived 
from subsets of the data. Excess observables also leave open the option of estimating additional parameters whose 
importance may become known only after the data are collected. For all these reasons, the experiment design usually 
provides for many more observables than the expected number of parameters. In fact, a typical astrometric/geodetic 
experiment will have several hundred observables per baseline. 

Given a network of stations, how arc the source observations chosen? To answer this question, and to generate 
an optimal observing schedule, the experimenter is guided by a few key principles. First the experiment must have 
sufficient spatial and temporal sampling in order to permit unambiguous separation of the parameters of interest 
during analysis. As an example of the wqqA for good spatial sampling, consider the dry zenith tropospheric delay 
parameter and the station local vertical coordinate v. The measured delays t will have sensitivities to these 
parameters that scale approximately as 



where E is the elevation angle above the horizon. These two parameters will be correlated if E is sampled only close 

to 45°. In order to separate the parameters the experiment must observe sources over a wide range of elevations. In 
particular the low-elevation {E < 10°) measurements will be most sensitive to the troposphere, while those at higher 
elevation angles will be most sensitive to the station vertical coordinate. 

Next we consider the need for good temporal sampling. Simply put, the experiment must sample the parameter 
of interest more quickly than it changes. For periodic effects this implies that one must take at least two samples 
within one full period of the effect of interest. In astrometric and geodetic VLBI one is often interested in tides with 
periods of 12 and 24 hours, 2 weeks, 1 month, a half a year, and one year. Components of the Earth's nutation have 
periods ranging from days up to 18.6 years. In practice experiments collect data at the rate of many observations per 
hour for a full 24-hour session. Such sessions have been carried out many times per year during the last two decades. 



dr/dZd oc 1/sini; 
dr/dv oc —sinE 



(2.1) 
(2.2) 
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allowing VLBI to measure parameters which vary on all of the aforementioned time scales. Day-long experiments also 
have the benefit of allowing observations through the full range of source right ascension. This minimizes problems 
that would otherwise arise in connecting observations from one right ascension zone to another. As a result, VLBI 
can avoid the zonal RA errors that have plagued optical observations which can only be performed at night. Because 
VLBI measurements are done at microwave frequencies, they are also much less sensitive to weather (clouds, rain) 
than optical techniques. 

The experimenter must next decide how long to observe each source. While collecting more data bits improves 
the detection signal-to- noise ratio (SNR), the desire to increase the SNR must be balanced against the reality of 
finite recording bandwidths which are 50-100 Mbits/sec for the Mark III VLBI system (Rogers et ai, 1983). As a 
result the optimal integration times typically range from 2 to 13 minutes. A technique known as bandwidth synthesis 
(Rogers, 1970) provides the highly accurate delays possible with a large radio frequency bandwidth, while requiring 
a recording bandwidth that is only a small fraction of the former. This technique has been essential in overcoming 
the limitations of currently available recording bandwidths. For very weak sources, much longer integration times 
may be required. The coherent integration time may be extended considerably with the phase referencing technique 
(Lestrade et ai, 1990), in which a strong nearby source is observed alternately with the weak source, and the phase 
of the strong source is used as a reference. 

Lastly, we consider how best to use a given network of antennas. For antennas separated by a few thousand 
kilometers or less it is often practical to have all antennas simultaneously observing the same source. With such modest 
baselines, each station can observe over a broad spatial extent. However, for networks with baselines significantly 
longer than an Earth radius, the desires for simultaneous visibility of the source and for broad geometric coverage 
start to be in conflict. Because the area of the sky that is simultaneously visible from both observing sites shrinks 
rapidly as baselines approach an Earth diameter in length, options for observing strategies become severely limited. 
The solution has been to create "sub-networks". This observing strategy assigns a subset of the stations in the 
network to observe one source, while another subset of stations observes a different source. While this complicates 
the experiment logistics, it has allowed the design of experiments with very strong geometries. 

Having considered the experiment design, we will now discuss how the signal is detected and recorded. For a 
Cassegrain antenna (see Fig. 9 in Sec. III.G), the incoming signal first strikes the primary paraboloidal dish of the 
antenna, is then reflected up to the hyperboloidal subreflector, and flnally into the feed horn on the central axis. For a 
prime focus antenna, the signal goes directly from the paraboloidal reflector to the feed at the focus position. Past the 
feed horn, the signal undergoes an initial stage of ampliflcation, either by transistor amplifiers such as field-effect or 
high-electron-mobility transistors or by a traveling wave maser. System temperatures typically range from 20-100 K 
at S and X band. After the first stage of amplification the signal is heterodyned from radio frequency to a lower 
intermediate frequency of several hundred MHz. Next the signal is further heterodyned down to "video" frequencies. 
There the signal is band limited to a width of a few MHz, sampled, and digitized. Lastly, it is formatted and recorded 
on digital video tape (see the lower half of Fig. 1). 

The signal is digitized at a minimum level of quantization to make the most effective use of the limited recording 
bandwidth. The digitized signals arc recorded on high capacity video tapes for a period of several minutes, yielding 
a total of wlO gigabits of data. In routine observing sessions lasting 24 hours, this procedure is repeated many times, 
with as many as a hundred distinct sources sampled at several distinct hour angles. Note that the independent station 
clocks miist remain well cnoiigh synchronized to obtain signal samples that will form a coherent interference pattern. 
The signal fiux density is on the order of 1 Jansky (Jy = 10^^^ W m^^Hz^^), necessitating antennas with large 
collecting areas and high sampling and recording rates, as well as highly sensitive and stable detectors and frequency 
standards. Pointing and mechanical characteristics of the antenna structure must be sufficiently responsive to permit 
movement between sources widely separated in the sky within a few minutes. 

The fundamental measurement in VLBI times the arrival of the wave front at the two ends of the baseline. In order 
to take full advantage of the timing precision possible with current frequency standards (e.g., hydrogen masers stable 
to 1 part in 10^^) care must be taken to calibrate phase shifts induced by the measurement instrumentation. Such 
phase shifts can corrupt the estimated phase and hence the group delay of the incoming signal. In order to correct for 
these instrumental phase shifts a technique known as phase calibration has been developed (Rogers, 1975; Fanselow, 
1976; Thomas, 1978; Sigman, 1987). This technique compensates for instrumental phase errors by generating a signal 
of known phase, injecting this signal into the front end of the VLBI signal path, and examining the signal's phase 
after it has traversed the instrumentation. This calibration signal is embedded in the broadband VLBI data stream as 
a set of low level monochromatic "tones" along with the signal from the radio source. Typically the calibration signal 
power is < 1% of the broadband VLBI signal level. The tones are extracted from the data at the time of subsequent 
processing. 
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B. Correlation: Producing the interference pattern 



After the signals at each antenna site are collected and recorded, the next step initiates the analysis of a VLBI 

experiment. The signals recorded at all participating antennas are combined pairwisc, thereby producing the interfer- 
ence pattern. This progression is shown in Fig. 1. The facilities for signal combination are called correlators. Several 
such installations are presently operating in the U. S. A. (Haystack Observatory, Jet Propulsion Laboratory, U. S. 
Naval Observatory, Very Long Baseline Array), Germany (Bonn), and Japan (Kashima). They consist of special- 
purpose signal processing hardware which is used to determine the difference in arrival times at the two stations by 
comparing the recorded bit streams. Moran (1976) and Thomas (1980; 1981; 1987) discuss the details of this process. 
In brief, the two bit streams representing the antenna voltages as functions of time t, Vi{t) and V2{t), are shifted in 
time relative to each other until their cross-correlation function R is maximized: 



1 

Rir) = ^ j dtV^{t)Vi{t-T) 



(2.3) 



T is the averaging interval, * denotes complex conjugate, and the time r corresponds to the difference in arrival times 
(the normalization is arbitrary for the purposes of our present discussion). 

Following Thomas (1987), a simplified diagram of the processing of data bits in a VLBI correlator is shown in 
Fig. 2. For purposes of this illustration we have only shown 8-bit averaging, one trial delay offset (known as a 
"lag"), artificially small delays, high fringe rates, and large amplitudes. Various fixed delays which occur in actual 
implementations are omitted in order to emphasize only the essential steps. 

Input bit streams 



Station 1 data 



station 2 data 



Delay 



Offset by model delays 
2 bits Delay = 1 bit 



^1 M^M-l-kM^ 



3-level cosine 



Bit— by— bit product 



3— level sine 



+ | - | + | | + | - | jj] 
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Cosine T, = +2 

bit-l 



Sine S = -2 

bit-l 



FIG. 2. Functional diagram of a VLBI correlator. The input bit streams are processed to the point of yielding the sine and 
cosine correlation sums. 



Starting at the top. Fig. 2 shows a short segment of the bit stream from each of two stations. The data were 
recorded using 1-bit sampling; when the signal voltage ^ > the bit is "+", and when V <Q the bit is "— ". Some 
of the bit locations have been left blank in this illustration, in order to show more clearly only the eight bits that 
are followed through all steps of the example. In the second step, the bit streams have been delayed by 2 bits for 
station 1 and by 1 bit for station 2. In reality these delays are typically on the order of 10^ bits, where one bit is 
recorded every 0.25 /is. In step 3, the two bit streams have been cross-correlated. The resulting product is "+" when 
the corresponding bits from the two streams match, and "— " when they differ. 

Since the Earth's rotational velocity induces Doppler shifts on the order of 10~^ times the observing frequency, 
VLBI observations at X band (8.4 GHz) will have their signals oscillating at several kHz. VLBI signals are typically 
very weak (only 1 bit in 10"^ or 10^ correlates between the two stations), and one must therefore average over many 
bits to detect the signal. To achieve this the signals must first be "counter-rotated" . This involves multiplying the 
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cross-correlated bits by sine and cosine waves to remove the above-mentioned Doppler shifts. Thus counter- rotation 
may also be thought of as a digital heterodyning of the cross-correlation signal from kHz to near zero frequency. In 
step 4, the bits are counter-rotated by multiplying them with 3-level approximations of sine and cosine waves. (This 
counter-rotation can also be performed prior to cross-correlating in step 3). Step 5: having been coimter-rotated, 
the signal may now be averaged over many bits in order to allow the weak VLBI signals to be detected. For clarity 
the diagram shows averaging over only eight bits. The resulting cosine and sine sums are then root-sum-squared and 
normalized to obtain the signal amplitude, and the arctangent (sine/cosine) is taken to yield the signal's phase. Note 
that for the purposes of this illustration the signal amplitudes are unrealistically large. 

Several aspects of the process described above distort the extracted signal and must be accounted for. First we 
note that 1-bit sampling causes a loss in SNR and introduces a 2/7r scale factor in the amplitudes (van Vleck and 
Middleton, 1966). The quantization of the delay to steps of one bit changes both the amplitude and phase, while 
the use of a 3-level approximation to the sine and cosine functions scales the amplitude. Thomas (1980; 1987) and 
Thompson, Moran and Swenson (1986) discuss these effects and their rectification in detail. 



C. Post correlation: Generating the observables 

The correlation process is carried out in parallel for many (typically 14) frequency channels, with each channel 
producing average amplitudes and phases every 1-2 seconds, as described in the previous section. These results are 
stored for later analyses with post-correlation software, which is discussed in detail by Lowe (1992). Briefly, such code 
performs many functions in order to prepare the data for the final processing by the modeling/estimation software. 
Some of these tasks are rejection of unreliable data by means of statistical tests, establishment of reference times to 
and frequencies loq for the observables, phase calibration via the tones that were injected soon after signal detection, 
and removal of any rudimentary model information that was introduced during correlation. The core task of the 
post-correlation software is to take the set of phase samples 0(wi, tj) from the various frequency channels iVi and times 
tj, and to fit the set of (j){LOi,tj) with three parameters: the phase 00, the group delay Tgd, and the phase rate fpd- To 
accomplish this, tlw^ sc^t of first Fourier transformed from the frequency and time domain to the delay 

and delay rate domain respectively. Figure 3 shows the result of this operation for a high-SNR observation. Next 
the delay and delay rate domain data are searched for peaks (at 763.784ns, — 0.3878 ps/s for the example in Fig. 3). 
This peak supplies the a priori values for a least squares fit which determines the observables 0o, Tgdi and Tpd (phase, 
group delay, and phase delay rate). 




FIG. 3. Fourier transform of the phase 4>(ui,tj) into the delay, delay rate domain. 

These phase-derived observables are determined (for phase (j) and circular frequency lu) from a bilinear least-squares 
fit to the measured phases ^{u),t): 



{0J,t) = (/)o{uJo,to) + -g^(UJ - UJo) + -g^{t - to) 



(2.4) 



where the phase, group delay, and phase rate are respectively defined as 
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Tpd = ?^0/w, 



Tgd 



duj ' 



(2.5) 



Formal uncertainties of each observable are also produced (for the example in Fig. 3, ar^^ = 0.010 ns, afp^ = 
0.0008ps/s). When data at two distinct frequencies are present {e.g., S and X bands), charged particle (iono- 
spheric) effects are removed by applying a simple model of dispersion, and combined (S/X) observables are formed 
(see Sec. V.A, Eq. 5.11). Problems with resolving 27r ambiguities in the signal phase over large distances usually 
preclude the direct use of the phase delay observable Tpj. This problem can be overcome in differenced observations 
of pairs of nearby sources, for which "phase connection" is possible (Marcaidc et al., 1994). Under routine conditions, 
the uncertainties in the delay observables arc on the order of 10 picoseconds. Numerous aspects of hardware, instru- 
mentation, and software contribute to this limit (Rogers, 1991). During correlation and post-correlation processing, 
the amount of data is compressed by a factor of 10^, from terabits to kilobits. 

To summarize, the interferometer is capable of producing four data types: phase delay, group delay, phase delay 
rate, and amplitude. The time rate of change of group delay cannot be measured accurately enough to be useful for 
geodetic or astrometric purposes. The delay models discussed in Sees. Ill, IV, and V are directly applicable to either 
group delay or to phase delay. Amplitudes are not usually modeled in astrometric/geodetic VLBI [for an exception 
see Chariot (1990a; 1990b)]. Phase delay rate models may be built from the delay models as discussed immediately 
below. 

The model for the time rate of change of phase delay must be constructed from either analytical or numerical time 
derivatives. If the latter route is taken, then only models of delay are needed. The model phase delay rate tpd(t) is 
approximated as a finite difference R by 



In the limit A — > 0, this expression for differenced phase delay approaches the instantaneous time rate of change of 

phase delay (fringe frequency) at time t. In practice, A must be large enough to avoid roundoff errors, but should 
also be small enough to allow i? to be a reasonably close approximation to the instantaneous delay rate. A suitable 
compromise appears to be a A in the vicinity of 0.1 second. Fortunately, the capability to model interferometer 
performance accurately is relatively insensitive to the choice of A over a fairly wide range of values. To be specific, if 
rpd(i) is expanded to third order in A, the numerically evaluated rate becomes 



Thus the error made in approximating the time derivative as R is proportional to fpdA^. This term amounts to a 
few times 10~^^ s/s for a 10,000 km baseline and A = 0.1 s. Comparisons of numerical and analytic derivatives have 
demonstrated equivalence of the two methods to better than 10~^^ s/s. 

D. Modeling and parameter estimation 

Prior to the detailed exposition in the next three sections of the theoretical modeling required by geodetic and 
astrometric VLBI, we present a brief summary of the VLBI modeling and parameter estimation process, its history, 
and key references. To form a perspective, wc start by listing in Table I the approximate relative importance of the 
various portions of the VLBI model. Each tabulated maximum delay is given in distance units, and represents the 
potential model error that woidd be caused by its complete omission from the theoretical model. Of the seven major 
components in this table, the purely geometric delay clearly dominates, since it can approach one Earth diameter. 
Aberration effects due to the Earth's orbital motion are four orders of magnitude smaller. The remaining components 
range over several orders of magnitude, with Earth orientation, antenna structure, instrumentation, and atmospheric 
effects being the most important. Most portions of the station coordinate and Earth orientation models become 
significant only in analyses of data extending over time spans much longer than a typical observing session of 24 
hours. The maximum delays arising from these models were therefore estimated for data extending over one year. 
Models of the last four major categories, on the other hand, can exhibit their total variation over the duration of one 
experiment, and their maximum contributions were therefore estimated over a time span of one day. 

We also attempt to provide an "error budget" for the VLBI model by tabulating rough estimates of modeling un- 
certainties in the last column of Table I. Extreme situations are considered for these estimates: Earth-sized baselines, 
extended sources, large antenna axis offsets, uncalibrated instrumentation, and observations at low elevations. It will 
be noted that many of the components can be modeled to accuracies that are several orders of magnitude better than 



R = [rpd(i + A) - rpd(i - A)]/(2A) . 



(2.6) 



R = fpd{t) + Tpd{t) A2 



/6. 



(2.7) 
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their maximum contributions. Thus, while this table indicates the general significance of the physical effects, their 
ranking does not necessarily reflect the uncertainties of their contributions to VLBI observables. Interactions between 
theory and experiment are quite important in interpretation of VLBI measurements: empirical analyses point out 
inadequate parts of the model, which can then be improved. Such iteration has taken place to varying degrees for 
different model components. Presently, each of the three basic aspects of a VLBI measurement (source, intervening 
medium, and receiver) are seen to impose limits to VLBI resolution and accuracy that are roughly of comparable 
magnitude. 

The effort expended in improvement of various model components has by no means been proportional to their 
relative contributions to the observables. For example, the ionosphere can be nearly totally eliminated as an error 
source if the measurements are performed at two widely separated frequencies. In recent years, much more attention 
has focused on station motion effects, despite their smaller size. While recent methods of treating Earth orientation 
and troposphere errors have produced substantial reductions of their contributions to the VLBI error budget, they 
nevertheless remain as two of the leading error sources in contemporary VLBI analyses. Two other important errors 
are due to incomplete modeling of the time-dependent structures of the radio sources and receiving antennas. Both 
may amount to tenths of a milliarcsecond (mas) (1 nrad, many mm). For day-long Earth-based measurements, the 
dominant error is due to the troposphere: commonly used models of atmospheric propagation can be incomplete at 
nearly the 10-20 mm level. In fact, present troposphere modeling is probably less accurate than the formal precision of 
the VLBI observables. Details of tidal motions on the Earth's surface have also not been fully characterized at the few 
mm level. For longer periods (on the order of 1 year), aperiodic time-dependent processes come into play. At present, 
these are likewise not well enough understood a priori at the nrad or mm level. The error budget of Table I leads 
to a (conservative) total modeling uncertainty of approximately 30 mm. More complete models of the Earth's tidal 
response and atmosphere, the physical processes in quasars, and the mechanical response of large antenna structures, 
will be required to realize the full potential accuracy of the VLBI technique and to reduce this limit to below 10 mm. 
It is hoped that this review provides a good point of departure for such model improvements. 

The a priori model is refined by estimating model parameter corrections which best fit the data. In astrometric 
and geodetic applications of radio interferometry, the values of the observables from many different radio sources 
are processed by a multiparameter least-squares estimation algorithm to extract refined model parameters. As the 
accuracy of the observables improves with advances in instrumentation, increasingly complete models of the delays 
and delay rates need to be developed, which will be the thrust of this review. 

Linear least-squares methods are employed to extract the best estimates of model parameters from the VLBI 
observables. Fortunately, nonlincaritics in this mathematical problem are small, and the a priori knowledge of 
most aspects of the model is sufficiently good to permit efficient estimation via linearized least squares. The only 
operational complication is caused by the sheer volume of data: more than two million pairs of delay and delay rate 
observations have been accumulated in experiments through 1997. This necessitates carehil bookkeeping procedures 
to keep computing requirements manageable when tens or hundreds of thousands of parameters are estimated in 
simultaneous analyses of all extant data. The most important of these is the separation of the parameters into "local" 
and "global" categories in order to minimize the need to invert very large matrices. For example, in processing data 
which extend over many observing sessions, the clock and atmosphere model parameters apply only to individual 
sessions, while station and soiirce coordinates might be determined by the entire data set. Specialized least-squares 
or filtering methods are usually employed in VLBI parameter estimation. Particular attention needs to bo paid to 
numerical stability, and to retaining correlation information in experiments that are widely separated in time. The 
"square root information filter" (SRIF) (Bierman, 1977) is one such formalism for achieving these goals. 

Severe limitations are imposed on the information content of the VLBI observables by the fact that the geometric 
delay measures only the scalar product of the baseline and source direction vectors. Thus the data determine the 
relative directions of the sources and baseline, and only the orientation of the whole Earth relative to the sources is 
strongly determined. Weak ties of the baseline to geographic features of the Earth are provided by effects such as 
the Earth's motion in its orbit (aberration), tidal displacements, and tropospheric delay. These fix the baseline in 
the terrestrial frame at an uncertainty level of a few tenths of an arcsecond. A more rigid link can only be provided 
by the introduction of a reference station with known coordinates {e.g., from satellite techniques) in the terrestrial 
reference frame, and the adoption of a conventional origin for celestial orientation {e.g., the right ascension of one 
reference source). 

A further complication in analysis is introduced by degeneracy in the parameter space of the multi-faceted model. 
Linear combinations of a subset of the parameters describing one portion of the model may produce identical changes 
in the observables to those described by parameters of another portion of the model. The tidal displacements of a 
global network of stations provide an illustrative example: certain linear combinations of station tidal motions are 
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equivalent to the three components of rotation of the Earth as a whole. All such potential degeneracies must be 
identified and accounted for in parameter estimation procedures. 

The account of the physical models presented here is intended to be helpful in providing a general understanding 
of the analysis of VLBI observations. Since modeling is intimately connected to the software that is used for data 
analysis, we take advantage of our own experience and adhere to the outline of the mathematical description of model 
implementation for the multiparameter estimation code "Masterfit/Modest" (Fanselow, 1983; Severs and Jacobs, 1996) 
in the succeeding sections. The theoretical foundation of this large-scale computer implementation was laid by J. G. 
Williams (1970a), who also developed many of the original algorithms. Software and model development has continued 
at the Jet Propulsion Laboratory since the 1970s. Independent code of similar scope and ancestry is the "Calc/Solve" 
package, which was developed at the Goddard Space Flight Center (GSFC) during the same time period (GSFC, 1981). 
This software has been used in the National Aeronautics and Space Administration (NASA) Crustal Dynamics and 
Space Geodesy Projects (GDP and SGP), National Oceanic and Atmospheric Administration (NCAA) International 
Radio Interferometric Surveying (IRIS), and National Earth Orientation Service (NEOS) analyses. It is an outgrowth 
of the original algorithms of Hinteregger et al. (1972) and Robertson (1975), and has been evolving predominantly 
at the Goddard Space Flight Genter. Numerous mutations of Galc/Solve over the years have resulted in at least 
three related, but not identical, packages used for VLBI analyses in Germany (Campbell, 1988), Japan (Kunimori 
et al, 1993), and Spain (Zarraoa et al, 1989). More recently, independent software has been developed in Ukraine 
(Yatskiv et al, 1991), Prance (Gontier, 1992), Norway (Andersen, 1995), and Russia (Petrov, 1995). There have been 
occasional comparisons between the various codes (Sovers and Ma, 1985; Gontier, 1992; Rius et al, 1992). These 
comparisons have been beneficial in exposing discrepancies, correcting software problems, and clarifying the physical 
models. Such discrepancies have been surprisingly minor, given the complexity of these large computer programs. 

III. GEOMETRIC DELAY MODELS 

The geometric delay is defined as the difference in time of arrival of a signal at two geometrically separate points 
which would be measured by perfect instrumentation, perfectly synchronized, if there were a perfect vacuum between 
the observed extragalactic or Solar-System source and the Earth-based instrumentation. For Earth-fixed baselines, 
this delay is at most the light time of one Earth radius (20 milliseconds) due to non-transparency of the Earth. It 
can change rapidly (by as much as 3 /is per second) as the Earth rotates. While VLBI experiments are occasionally 
carried out with more than ten participating stations, the correlator generates observable delays and their time rates 
of change independently for each baseline connecting a pair of stations. Without loss of generality, the delay model 
can thus be developed for a single baseline involving only two stations, which is what we present here. The geometric 
component is by far the largest component of the observed delay. The main complexity of this portion of the model 
arises from the numerous coordinate transformations that are necessary to relate the celestial reference frame used 
for locating the radio sources to the terrestrial reference frame in which the station locations are represented. 

We use the term "celestial reference frame" to denote a reference frame in which there is no net proper motion 
of the extragalactic radio sources which are observed by the interferometer. This is only an approximation to some 
truly "inertial" frame. Currently, this celestial frame is a Solar-System-baryccntric (SSB), equatorial frame with the 
equator and equinox of Julian date 2000 January 1.5 (J2000.0) as defined by the 1976 International Astronomical Union 
(lAU) conventions, including the 1980 nutation series (Seidelmann et al, 1992; Kaplan, 1981). In this equatorial frame, 
some definition of the origin of right ascension must be made. The right ascension is nearly arbitrary (neglecting 
small aberrational, tidal, and gravitational bending effects due to the Earth's changing position within the Solar 
System). Any definition differs by a simple rotation from any other definition. The important point is that consistent 
conventions must be used throughout the model development. The need for this consistency has recently led to a 
definition of the origin of right ascension in terms of the positions of extragalactic radio sources (Ma et al, 1997). 
Interferometric observations of both natural radio sources and spacecraft at planetary encounters will then connect 
the optical reference frame and planetary ephemeris with the radio reference frame (Newhall et al, 1986; Folkner et 
al, 1994; Lestrade et al, 1995). 

Unless otherwise stated, we will mean by "terrestrial reference frame" a frame tied to the mean surface features of the 
Earth. The most common such frame, which we use here, is a right-handed version of the Conventional International 
Origin (CIO) reference system with the pole defined by the 1903.0 pole. In practice, the tie is most simply realized by 
defining the position of one of the interferometric observing stations, and then determining the positions of the other 
stations under this constraint. This tic further requires that the determinations of Earth orientation agree on the 
average with measurements of the Earth's orientation by the International Earth Rotation Service (lERS) (1996a) 
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[and its predecessor, Bureau International de I'Heure (BIH) (1983)] over some substantial time interval (« years). Such 
a procedure, or its functional equivalent, is necessary to tie the measurement to the Earth, since the interferometer is 
sensitive only to the baseline vector. With the exception of minor tidal and tropospheric effects, the VLBI technique 
does not have any preferred origin relative to the structure of the Earth. The rotation of the Earth does, however, 
provide a preferred direction in space which can be associated indirectly with its surface features. 

In contrast, geodetic techniques which involve the use of artificial satellites or the Moon are sensitive to the center 
of mass of the Earth (e.g. Vigue et at, 1992) as well as its spin axis. Thus, such techniques require only a definition of 
the origin of longitude. Laser ranging to the retroreflectors on the Moon allows a realizable definition of a terrestrial 
frame, accurately positioned relative to a celestial frame which is tied to the planetary ephemerides (Folkner et 
ai, 1994). The required collocation of the laser and VLBI stations is being provided by Global Positioning Satellite 
(GPS) measurements of baselines between VLBI and laser sites starting in the late 1980s {e.g. Ray et al., 1991). 
Careful definitions and experiments of this sort are required to realize a coordinate system of centimeter accuracy. 

Except for sub-centimeter rclativistic complications caused by the locally varying Earth potential (as discussed 
below), construction of the VLBI model for the observed delay can be summarized in 7 steps as: 

1. Specify the proper locations of the two stations as measured in an Earth- fixed frame at the time that 
the wave front intersects station 1. Let this time be the proper time t[ as measured by a clock in the 
Earth-fixed frame. 

2. Modify the station locations for Earth-fixed effects such as solid Earth tides, tectonic motion, and other 
local station and antenna motion. 

3. Transform these proper station locations to a geocentric celestial coordinate system with its origin at the 
center of the Earth, and moving with the Earth's center of mass (but not rotating). This is a composite 
of 12 separate rotations. 

4. Perform a Lorentz transformation of these proper station locations from the geocentric celestial frame to a 
frame at rest relative to the center of mass of the Solar System, and rotationally aligned with the celestial 

geocentric frame. 

5. In this Solar-System-barycentric (SSB) frame, compute the proper time delay for the passage of the 
specified wave front from station 1 to station 2. Correct for source structure. Add the effective change in 
proper delay caused by the differential gravitational retardation of the signal within the Solar System. 

6. Perform a Lorentz transformation of this SSB geometric delay back to the celestial geocentric frame moving 
with the Earth. This produces the adopted model for the geometric portion of the observed delay. 

7. To this geometric delay, add the contributions due to clock offsets (Sec. IV), to tropospheric delays, and 
to the effects of the ionosphere on the signal (Sec. V). 

As indicated in step 5, the initial calculation of delay is carried out in a frame at rest relative to the center of mass 
of the Solar System (SSB frame). First, however, steps 1 through 4 are carried out in order to relate proper locations 
in the Earth-fixed frame to corresponding proper locations in the SSB frame. Step 4 in this process transforms station 
locations from the geocentric celestial frame to the SSB frame. This step incorporates spccial-relativistic effects to all 
orders of the velocity ratio v/c. From a general relativistic point of view, this transformation is a special relativistic 
transformation between proper coordinates of two local frames (geocentric and SSB) in relative motion. For both 
frames, the underlying gravitational potential can be taken approximately as the sum of locally constant potentials 
caused by all masses in the Solar System. The complications caused by small local variations in the Earth's potential 
are discussed below. The initial proper delay is then computed (step 5) in the SSB frame on the basis of these SSB 
station locations and an a priori SSB source location. A small proper-delay correction is then applied to accoimt 
for the differential gravitational retardation introduced along the two ray paths through the Solar System, including 
retardation by the Earth's gravity. A final Lorentz transformation (including all orders of v/c) then transforms the 
corrected SSB proper delay to a model for the observed delay in the celestial geocentric reference frame. Note that, 
rather than transforming between frames, some formulations account for special relativistic effects as an ad hoc term 
expanded to order (w/c)^. While this may be convenient for calculating the delay, it obscures the important physical 
insights that are gained by a clear definition of the reference frame used for each step. 

Since the Earth's gravitational potential Ue varies slightly through the Earth {AUe/c^ « 3.5 x 10^^° from center 
to surface), the specification of proper distance for a baseline passing through the Earth is not as straightforward 
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with respect to the Earth's potential as it is with respect to the essentially constant potentials of distant masses. To 
overcome this difficulty, VLBI-derived station locations are now customarily specified in terms of the "TDT (Terrestrial 
Dynamic Time) spatial coordinates" that are used in Earth-orbiter models. A proper length that corresponds to a 
modeled baseline can be obtained through appropriate integration of the local metric (Shahid-Saless et al., 1991). 
Such proper lengths deviate slightly (< 3 mm) from baselines modeled on the basis of the TDT convention in the 
worst case (a full Earth diameter). In practice, such a conversion is not necessary if baseline measurements obtained 
by different investigators are reported and compared in terms of TDT spatial coordinates. 

The model presented here has been compared (Treuhaft and Thomas, 1991) with the "1-picosecond" relativistic 
model for VLBI delays developed by Shahid-Saless et al. (1991). When reduced to the same form, these two models 
are identical at the picosecond level, term by term, with one exception. (Treuhaft and Thomas, 1991) show that a 
correction is needed to the Shahid-Saless et al. modeling of the atmospheric delay in the SSB frame. This correction 
changes the Shahid-Saless et al. result by as much as 10 picoseconds. Comparisons of this formulation with the 
Goddard Space Flight Center implementation of the "consensus" model (Eubanks, 1991; lERS, 1996b) show agreement 
to better than 1 ps, suggesting that the atmospheric delay problem has been solved in the consensus model. The 
remainder of this section provides the details for the first six steps of the general outline given above. 

A. Wave front arrival time dIfTerence 

The fundamental part of the geometric model is the calculation (step 5 above) of the time interval for the passage of 
a wave front from station 1 to station 2. This calculation is actually performed in a coordinate frame at rest relative 
to the center of mass of the Solar System. The SSB coordinate system is used because source positions, planetary 
orbits, and motions of interplanetary spacecraft are most naturally modeled in this frame. This part of the model is 
presented first to provide a context for the subsequent sections, all of which are heavily involved with the details of 
time definitions and coordinate transformations. We will use the same subscript and superscript notation which is 
used in Sec. III.E to refer to the station locations as seen by an observer at rest relative to the center of mass of the 
Solar System. 

First, we calculate the proper time delay that would be observed if the wave front were planar. This calculation 
is next generalized to a curved wave front, and finally we take into account the incremental effects which result from 
the wave front propagating through the various gravitational potential wells in the Solar System. 

1. Plane wave front 

Consider the case of a plane wave moving in the direction k with station 2 having a velocity ^2 shown in Fig. 4. 
As mentioned above, distance and time are to be represented as proper coordinates in the SSB frame. The speed of 
light c is set equal to 1 in the following formulation. The proper time delay is the time it takes the wave front to move 
the distance I at speed c. This distance is the sum of the two solid lines perpendicular to the wave front in Fig. 4: 



Note that station 2 has moved since ti : the second term represents the distance that station 2 moves before receiving 
the signal at t2. This leads to the following expression for the geometric delay: 



The baseline vector in the SSB frame, r2{ti) — ri(ti), is computed from proper station locations Lorentz transformed 
from the geocentric celestial frame, using Eq. (3.158) in Sec. III.E. The source unit vector in the direction of signal 
propagation (source to receiver) 



t2-ti=k- [r2(ti) - ri(ti)] + k • /32[t2 - ti 



(3.1) 



t2-tl = 



k- [r2(ti)~ri(ti)] 

1 - k • /32 



(3.2) 




(3.3) 



is calculated from the SSB frame angular source coordinates a (right ascension) and 6 (declination). 
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POSITION OF STATION #2 
WHEN WAVE FRONT CROSSES 
IT AT TIME 



POSITION OF STATION #1 
WHEN WAVE FRONT CROSSES 
IT AT TIME 



FIG. 4. Geometry for calculating the transit time of a plane wave front. Some refinements have been added to the schematic 
diagram of Fig. 1. 



2. Curved wave front 

In the case of a signal generated by a radio source within the Solar System it is necessary to include the ef- 
fect of the curvature of the wave front. As depicted in Fig. 5, let a source irradiate two Earth- fixed stations 
whose positions are given by ri,2{t) relative to the Earth's center. The position of the Earth's center in the 
SSB frame (Rc(ti)) as a function of signal reception time ti at station 1 is measured relative to the position of 
the emitter at the time of emission (te)- While this calculation is actually done in the Solar System barycen- 
tric coordinate system, the development that follows is by no means restricted in applicability to that frame. 



FIG. 5. Geometry for calculating the transit time of a curved wave front. It is only necessary to consider curvature effects for 
sources closer than «30 light years. 

Suppose that a wave front emitted by the source at time tg reaches station 1 at time ti and arrives at station 2 at 
time t2- The geometric delay in this frame will be given by 




18 



T = i2-ii = |R2(t2)|-|Rl(tl)| , (3.4) 

where all distances are again measured in units of light travel time. If we approximate the velocity of station 2 by 



and use the relation (i=l,2) 



we obtain: 



^ t2-ti ^ ' 



Ri(ti) = Rc(ti) + ri(fi) (3.6) 



r = |Rc(ii) + r2(ii) + f3^T\ - |R,(ii) + ri{ti)\ 
= Rc{h) [ |Rc + £2| - |Rc + £i| ] , (3.7) 



where 



£i=ri{ti)/Rc{ti) 

£2 = [r2(ii)+/92^]/i2c(ti) • (3.8) 

For £i and £2 < 10~^, only terms of order £^ need to be retained in expanding the expression for r in Eq. (3.7) for 
computational purposes on machines that employ 15-16 decimal digit representations. This gives: 

^ ^ Rc- [r2(^i)-ri(^i)] ^ i?cAe(r) 

[1 - Re-/32] 2 [1 - Re-/32] ' 

where to order 

Ac(r) = [el - ej] - [(Rc-£2)' + (Rc-£i)' + (Rc-£2)^ - (Rc-£2)£2 - (Rc-£i)^ + (Rc-£i)£?] • (3.10) 

The first term in Eq. (3.9) is just the plane wave approximation, i.e., as Rc — » 00, Rc k, with the second term in 
brackets in Eq. (3.10) approaching zero as /Rc- Given that the ratio of the first term to the second term is « r/Rc, 
wave front curvature is not calculable using sixtccn-dccimal-digit arithmetic if i? > lO^^r. For Earth-fixed baselines 
that are as long as an Earth diameter, requiring that the effects of curvature be less than 1 ps » 0.3 mm implies that 
the above formulation (Eq. 3.10) must be used for i? < 3 x 10^^ km, or approximately 30 light years. At the same 
accuracy level, fourth and higher order terms in e become important for R < 2 x 10'' km, or approximately inside the 
Moon's orbit. Fukushima (1994) has presented a formulation which should be applicable for Earth-based VLBI with 
radio sources as close as the Moon. 

The procedure for the solution of Eq. (3.9) is iterative for e < 10~^, using the following: 

ficAc(T„-i) 

= "° + 2[1 - R./32] ' ^'-''^ 

where 

'^0 — '^plane wave • (3.12) 

For £ > 10^^, directly iterate on Eq. (3.7) itself, using the procedure: 

Tn = i?c|Rc + e2(r„-i)| - i?c|Rc + £i| , (3.13) 
where again to is the plane wave approximation. 
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3. Gravitational delay 



As predicted by Einstein's theory of general relativity (Einstein, 1911, 1916), an electromagnetic signal propagating 

in a gravitational potential is retarded relative to its travel time in field-free space. Shapiro (1964; 1967) pointed 
out that this produces both a deviation from a straight-line path and a time delay. For VLBI this implies that the 
computed differential arrival time value of the signals at ri(ti) and r2(t2) must be corrected for gravitational effects. 
For Earth-based experiments the dominant contribution comes from the Sun. Gravitational potential effects and 
curved wave front effects can be calculated independently of each other since the former are a small perturbation 
(pa 8.5 microradians or < 1 .75), even for Sun-grazing rays. 

Our formulation of the relativistic light travel time is based on Moyer's (1971) implementation of the work of 
Tausner (1966) and Holdridge (1967). For the (exaggerated) geometry illustrated in Fig. 6, the required correction to 
coordinate time delay due to the pth gravitating body is 



A (1 +7ppn)Mp 

= :3 



Ts + r2(^2) +rs2 \ ''s +ri{ti) + rsi 



(3.14) 



Vs + r'2ih) - J \rs + ri{ti) - rsi J 
where Vgi {i = 1, 2) is the separation of station i from source s: 

rsi = \ri{U) - r,{te)\ . (3.15) 

The gravitational constant Hp is 

Hp = Gnip , (3.16) 

where G is the universal gravitational constant, and rrip is the mass of body p. A higher-order term (oc iip^ /c") 
contributes to Aqp only for observations extremely close to the Sun's limb (lERS, 1992). Here 7ppN is the 7 factor 
in the parametrized post-Newtonian gravitational theory [e.g. Misner et al., 1973). In general relativity 7pp^ = 1, 
but it can be allowed to be an estimated parameter to permit experimental tests of general relativity, as originally 
proposed by Shapiro (1964). A dramatic demonstration of the importance of general relativity in estimating the 
delay can be made by "turning off" this added delay by setting 7pp^ = — 1 in the model. This can increase the 
discrepancy between theoretical and observed delays by an order of magnitude because the Sun's gravity induces 
a bending of 20 nrad (4 milliarcseconds) even for ray paths 90° removed from its position relative to the Earth. 




FIG. 6. Schematic representation of the geodesic connecting the source and receiver in the presence of a gravitational mass. 



Depending on the particular source- receiver geometry in a VLBI experiment, a mimber of useful approximations 
are possible for the correction Aqp of Eq. (3.14). Dropping the time arguments in Eq. (3.14), we have: 
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A _ (-*- + 7ppn)Mp 1 

^Gp = 3 in 



rs + r2+ rs2 \ I + ri - r^i 
rs + ri+ r^i / 1 + ra - 



(3.17) 



This formulation is appropriate for the most general geometry, in which « « r^j. For the practical case of 
Earth-based VLBI with distant sources and closely spaced VLBI receivers, however, |r2 — ri|/ri 0,rj/rs 0. The 
gravitational time delay Acp may then be expanded in terms of Vi/rg, Vsi/vs- Making use of the relationship 



r,i = [rl - 2r, ■ + r?]'/" « - • r, (3.18) 



leads to 



for ri/vs — > 0. If we define r2 = ri + Ar, and require that |Ar|/ri — > (short baselines) then 

r2(l + r2-?.)wri(l+?i-?,) + Ar •(?!+?,) . (3.20) 
Substituting into Eq. (3.19) and expanding the logarithm, we obtain: 

A -(1 + 7ppn)Mp (''2 - ri) ■ (ri +r,) 

Agp = 5 7T~r^ — =^1 • (<5-^l) 

ri(l + ri-rs) 

These three formulations of Eqs. (3.17, 3.19 or 3.21) are computationally appropriate for different situations. For 
Earth-based VLBI observations, the correction Acp is calculated for each of the major bodies in the Solar System 
(Sun, planets, Earth, and Moon). The immense central mass of the Galaxy (on the order of 10^^ solar masses) 
contributes an additional delay. Assuming a mass concentration of 2/3 in the nucleus and the current estimate of 
our distance from the center of 8.5 kpc (Binney and Tremaine, 1987), it can be estimated that the gravitational 
influence of the Galactic center on electromagnetic signals exceeds that of the Sun by a factor of «40. Geometrically, 
this causes a bending of 4 arc seconds when ray paths approach within « 10° of the center (for the closest routinely 
observed radio source). However, because of the very slow rotation of the Solar System about the Galactic center 
(240 million years), the time variation of this bending is extremely small. Therefore, for practical purposes. Galactic 
bending merely produces a quasi-static distortion of the sky, and can be absorbed into the source coordinates which 
form the basis of our model of the celestial sphere. 

Before the correction Aqp can be applied to a proper delay computed according to Eq. (3.2), it must be converted 
from a coordinate-delay correction to a proper-delay correction appropriate to a near-Earth frame. For such proper 
delays, the gravitational correction is given {e.g. Hellings, 1986) to good approximation by 

A^p = AGp-(l+7ppJt/T, (3.22) 

where r is the proper delay given by Eq. (3.2), and where U = p^p/r.pC? is the negative of the gravitational potential 
of the given mass divided by c^, as observed in the vicinity of the Earth {U is a positive quantity). The Ut term is 
a consequence of the relationship of coordinate time to proper time, and the 7ppNJ/r term is a consequence of the 
relationship of coordinate distance to proper distance. 
The total gravitational correction used is 

JV 

A^ = E^Gp. (3.23) 

p=i 

where the summation over p is for the major bodies in the Solar System. For the Earth, the (1 + ')-p-p^)UT term in 
Eq. (3.22) is omitted if one wishes to work in the "TDT spatial coordinates" that are used in reduction of Earth-orbiter 
data. The scale factor (1 + ^ppf^)U is approximately 19.7 x 10~® for the Sun. A number of other conventions are 
possible. One of these, which includes the (1 + 7pp^- )C/t term for the Earth evaluated at the Earth's surface, yields 
an additional scale factor of 1.4 x 10~^. In either case, the model delay is decreased relative to the proper delay. 
Consequently, all inferred "measured" lengths increase by the same fraction relative to proper lengths [i.e., by 19.7 
parts per billion and 21.1 ppb in the two cases). 
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Some care must be taken in defining the positions given by r^, r2{t2), and ri(ii). We have chosen the origin to 
be the position of the gravitational mass at the time of closest approach of the received signal to that object. The 
position Fg of the source relative to this origin is the position of that source at the time te of the emission of the 
received signal. Likewise, the position ri{ti) of the ith receiver is its position in this coordinate system at the time of 
reception of the signal. Even with this care in the definition of the relative positions, we are making an approximation, 
and implicitly assuming that such an approximation is no worse than the approximations used by Moyer (1971) to 
obtain Eq. (3.14). 

Some considerations follow, concerning the use of appropriate times to obtain the positions of the emitter, the 
gravitational object, and the receivers. For a grazing ray emitted by a source at infinity, using the position of the pth 
gravitating body at the time of reception of the signal at station 1 (i,.) rather than at the time of closest approach 
of the signal to mp (ta) can cause a substantial error on baselines with dimensions of the Earth, as shown by the 
following calculation. From Fig. 7, the distance of closest approach R changes during the light transit time tti of a 
signal from a gravitational object at a distance Reg by 

AR « Reg a itr = aRlc/c , (3.24) 

where a is the time rate of change of the angular position of the deflecting body as observed from Earth. Since the 
deflection is 

6^2 ^^ +^7'^^^^ I ^ I (3.25) 



it changes by an amount 



during the light transit time. 




Ae = -e|^)=-e(^) (3.26) 




FIG. 7. Schematic representation of the motion of a gravitating object (Sun) during the transit time of a signal from the point 
of closest approach to reception by an antenna on Earth. Both move between the time ta at closest approach R{ta) and the 
time tr at reception. This motion needs to be taken into account for the highest accuracy. 

Wc consider the two bodies of largest mass in the Solar System: the Sim and Jupiter. For rays that just graze 
the surfaces, their respective deflections AQ are 8470 and 79 nanoradians. From Solar System ephemerides, the 
barycentric angular velocities a are estimated to be w 0.05 and 17 nrad/s for the Sun and Jupiter. (The Sun's motion 
in the barycentric frame has approximately the orbital period of Jupiter, ft!l2 years, with a radius on the order of the 
Sun's radius). Using approximate radii and distances from Earth to estimate Reg and O, Eq. (3.26) gives 30 nrad 
for Jupiter; the corresponding value for the Sun is only 0.05 nrad. For a baseline whose length equals the radius of 
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the Earth, 5(A9)i?E is thus approximately 200 and 0.3 mm for Jupiter and the Sun, respectively. The effect is much 
smaller for the Sun in spite of its much larger mass, due to its extremely slow motion in the barycentric frame. 

In view of the very rapid decrease of gravitational deflection with increasing distance of closest approach, it is 
extremely improbable that a random VLBI observation will involve rays passing close enough to a gravitating body 
for this correction to be of importance. Exceptions are experiments that are specifically designed to measure planetary 
gravitational delay (e.g. Treuhaft and Lowe, 1991). In order to guard against such an unlikely situation in routine 
work, and to provide analysis capability for special experiments, it is prudent to perform the transit-time correction for 
all planets for all observations. To obtain the positions of the gravitational objects, we employ an iterative procedure, 
using the positions and velocities of the objects at signal reception time. If R(fr) is the position of the gravitational 
object at signal reception time tr then that object's position R(ta) at the time ta of closest approach of the ray path 
to the object Ra was 

IL{ta) = R{tr) - {tr - ta)V , (3.27) 
tr-ta = |Rc|/c . (3.28) 

This correction is done iteratively, using the velocity Y{tr) as an approximation of the mean velocity V. Because 
u/c « 10~^, the iterative solution, 

R„(ta)=R(tr) - |Rn-l(ta)| V(i,)/C (3.29) 

rapidly converges to the required accuracy. 



B. Time information 

Before continuing with the description of the geometric model, some definitions must be introduced concerning time- 
tag information in the experiments, and the time units which will appear as arguments below. A general reference 
for time definitions is the Explanatory Supplement (Seidelmann et a/., 1992). The data epoch for VLBI observables 
is taken from the UTC (Coordinated Universal Time) time tags in the data stream at station 1, UTCi. This time is 
converted to Terrestrial Dynamic Time TDT, with the conversion consisting of the following components: 

TDT (TDT - TAI) + (TAX - UTCiers) + (UTCiers - UTCstd) 

+ (UTCsTD - UTCi) + UTCi . (3.30) 

The four offsets in Eq. (3.30) thus serve to convert the station 1 time tags to TDT. Their meanings are as follows: 

1. TDT— TAI is 32.184 seconds by definition; TAI (Temps Atomique International) is atomic time. 

2. TAI— UTCiers is the offset between atomic and coordinated time. It is a published integer second offset 
(accumulated leap seconds) for any epoch after 1 January, 1972. Prior to that time, it is a more complicated 
function, which will not be discussed here since normally no observations previous to the mid-1970s 
are modeled. The International Earth Rotation Service (lERS), its predecessor, Bureau International 
de I'Heure (BIH), and Bureau International des Poids et Mesures (BIPM) are the coordinating bodies 
responsible for upkeep and publication of standard time and Earth rotation quantities. 

3. UTCiers— UTCstd is the offset in UTC between the coordinated time scales maintained by the lERS 
(BIPM) and secondary standards maintained by numerous national organizations. For VLBI stations in 
the U. S. these secondary standards are those of the National Institute of Standards and Technology in 
Boulder, Colorado, and the U. S. Naval Observatory in Washington, DC. These offsets can be obtained 
from BIPM Circular T {e.g. BIPM, 1997). 

4. UTCstd^ UTCi is the (often unknown) offset between UTC kept by station 1 and the secondary national 
standard. This may be as large as several microseconds, and is not precisely known for all experiments. 
It can be a source of modeling error: an error At in epoch time causes an error of w BuiEAt = 7.3 x 10~^ 
mm per km of baseline B per /is of clock error At, where cje is the rotation rate of the Earth. Even in 
the extreme case of a 10,000 km baseline, however, this amounts to only 0.7 mm per /xs, while present-day 
clock synchronization is usually at least an order of magnitude better, at the 0.1-/is level. 
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Coordinated Universal Time and TDT are not convenient for modeling Earth orientation. This is normally done in 
terms of UTl, a quantity that is proportional to the angle of rotation of the Earth. The a priori offset UTl— UTC and 
the position of the Earth's rotation pole can be obtained by interpolation of the lERS Bulletin A smoothed values. 
However, any other source of UTl— UTC and pole position could be used provided it is expressed in a left-handed 
coordinate system (see Sec. III.D.l). Part of the documentation for any particular set of results needs to include a 
clear statement of what values of UTl— UTC and pole position were used in the data reduction process. 

For the Earth model based on the lAU conventions, the following definitions are employed throughout (Kaplan, 
1981): 

1. Juhan date at epoch J2000.0 = 2451545.0. 

2. All time arguments denoted by T below are measured in Julian centuries of 36525 days of the appropriate 
time relative to the epoch J2000.0, i.e., T = ( JD - 2451545.0)/36525. 

3. For the time arguments used to obtain precession and nutation, or to refer to the Solar System ephemeris, 
Barycentric Dynamic Time (TDB, Temps Dynamique Barycentrique) is used. This is related to Terrestrial 
Dynamic Time (TDT, Temps Dynamique Terrestre) by the following approximation: 

TDB = TDT + 0.^001658 sm{g + 0.0167 sin(£f)) , (3.31) 

where 

g = 27r(357.°528 + 35999.°050 r)/360° (3.32) 

is the mean anomaly of the Earth in its orbit. 

For present analyses of Earth-based VLBI observations, Eq. (3.31) is adequate. Moyer (1981) gave a more accurate 
relation between TDB and TDT by accounting for the dominant gravitational effects of the major planets. Hirayama 
et al. (1987) and Fairhead and Bretagnon (1990) have extended this work to the nanosecond level. In the future, 
TDT and TDB will be replaced by two new time scales, TCG and TCB, geocentric and barycentric coordinate time 
(Pukushima at al, 1986). These will eliminate the need to rescale spatial coordinates that was discussed in Sec. III.A.3. 



C. Station locations 

Coordinates of the observing stations are expressed in the Conventional International Origin (CIO 1903) reference 
system, with the reference point for each antenna defined as in Sec. III.G. The present accuracy level of space geodesy 
makes it imperative to account for various types of crustal motions. Among these deformations are tectonic motions, 
solid Earth tides, ocean effects, and alterations of the Earth's surface due to local geological, hydrological, and 
atmospheric processes. Mismodeled effects will manifest themselves as temporal changes of the Earth-fixed baseline. 
It is therefore important to model all crustal motions as completely as possible. The current level of mismodeling of 
these motions is probably one of the leading sources of systematic error (along with the troposphere) in analyses of 
VLBI data. 

Evaluation of the time dependence of station locations is most simply done by estimating a new set of coordinates in 
the least-squares process for each VLBI observing session {e.g., Fig. 19 in Sec. VLB). Subsequent fits to these results 
can then produce estimates of the linear time rate of change of the station location. For rigorous interpretation of 
the statistical significance of the results, correlations of coordinates estimated at different epochs must be fully taken 
into account when calculating uncertainties of the rates. The advantage of this approach is that the contribution 
of each session to the overall time rate may be independently evaluated, since it is clearly isolated. Any nonlinear 
{e.g., seasonal) variations are easily detected. Also, no tectonic model information is imposed on the solution. An 
alternative approach is to model long-period tectonic motion directly, by introducing time rates of change of the 
station coordinates as parameters. The model is linear, with the position of station i at time t, rj = {xi,yi,Zi), 
expressed in terms of its velocity as 

ri = r°+ri(i-io) • (3.33) 
Here to is a reference epoch, at which the station position is r? = (x?, y^,z^). 
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1. Tectonic plate motion 



As alternatives to estimating linear time dependence of the station coordinates from VLBI experiments, several 

standard models of tectonic plate motion arc available. They all describe the motion as a rotation of a given rigid 
plate (spherical cap) about its rotation pole on the surface of a spherical Earth. Time dependence of the Cartesian 
station coordinates of station i which resides on plate j is expressed as 

ri = r° + (u;^-xrO)(i-io) , (3.34) 

where u}^ = {u^, co^, w^) is the angular velocity vector of the plate. 

These tectonic motion models are based on paleomagnetic data spanning millions of years, but they also provide an 
excellent quantitative characterization of present-day plate motions. This attests to the smooth character of tectonic 
motion over immense time scales. Although the possibility of slow relative motions of the Earth's land masses was first 
suggested by Wegener before World War I (Wegener, 1912) mantle convection was identified as the driving mechanism 
only much later, and the reality of this phenomenon was not generally accepted until the 1960s (Menard, 1986). The 
first quantitative global tectonic motion model is due to Minster and Jordan (1978), and was also the first to be used 
in VLBI analyses. It is denoted AMO-2 in the original paper. More recent models, denoted Nuvel-1 and NNR-Nuvel-1, 
are due to DeMets et al. (1990) and Argus and Gordon (1991), respectively. In Nuvel-1, the Pacific plate is stationary. 




FIG. 8. Schematic directions and magnitudes of tlic motions of the major tectonic plates in the NNR-Nuvel-IA model. The 
longest arrows represent motions of approximately 10 cm/yr. 

while NNR-Nuvel-1 is based on the imposition of a no- net-rotation (NNR) condition. Inadequate knowledge of the 

internal mechanics of the Earth makes uncertain any absolute determination of plate rotation relative to the deep 
interior. Hence the NNR condition is customarily imposed: v x r integrated over the Earth's surface is constrained to 
be zero, where v is the velocity at point r on one of the rigidly rotating tectonic plates. With some notable exceptions, 
the Nuvel-1 models give rates that arc very close to those of the AMO-2 model. The AMO-2 India plate has been split 
into two: Australia and India, and there are four additional plates: Juan de Fuca, Philippine, Rivera and Scotia. A 
recent revision of the paleomagnetic time scale has led to a rescaling of the Nuvel-1 rates. These "Nuvel-IA" and 
"NNR-Nuvel-1 A" model rates are equal to the Nuvel-1 and NNR-Nuvel-1 rates, respectively, multiplied by a factor 
of 0.9562 (DeMets et al., 1994). Table II shows the angular velocities of the 16 tectonic plates in the NNR-Nuvel-IA 
model. Thirteen major plate motions are also shown in Fig. 8, with the velocities depicted approximately to scale. 
The largest arrows represent motions of « 10 cm/yr. A more detailed example of the effect of tectonic motion on 
relative station locations is given in Fig. 19 of Sec. VLB, which shows the convergence of the North American and 
Pacific plates. 
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2. Tidal station motion 

The Earth is not perfectly rigid, and its crust deforms in response to the gravitational attraction of other massive 

bodies. The most important of these are the Sun and Moon, whose periodic orbits induce time- varying displacements 
of the Earth's crust. Such motions, with periodicities ranging from hours to years, are called tidal effects. Tides 
produce station displacements that can be far larger than those caused by tectonic motions in a year, necessitating 
their inclusion in models of VLBI observables. The tidal displacements can be classified into several categories, of 
which contemporary VLBI models normally include four. In the standard terrestrial coordinate system, these tidal 
effects modify the station location ro by an amount 

A = Asol + Apol + Aocn + Aatm , (3.35) 

where the four terms are due to solid Earth tides (the direct primary effect), pole tide, ocean loading, and atmosphere 
loading (secondary effects), respectively. Other Earth-fixed effects {e.g., glacial loading) can be incorporated by 
extending the definition of A to include additional terms. All four tidal effects arc most easily calculated in some 
variant of a VEN (Vertical, East, North) local geodetic coordinate system. Application of the transformation VW 
(given in the next section, Eqs. 3.61) transforms them to the geocentric Earth-fixed coordinate frame. 



a. Solid Earth tides 

Calculating the shifts of the positions of the observing stations caused by solid Earth tides is rather complicated 
due to the solid tides' coupling with the ocean tides, and the effects of local geology. Some of these complications 
are addressed below {e.g., ocean loading). An isolated simple model of Earth tides is the multipole response model 
developed by Williams (1970b), who used Melchior (1966) as a reference. Let Rp be the position of a tide-producing 
body, and ro the station position. To allow for a phase shift ip of the tidal effects from a nominal value of 0, the phase- 
shifted station vector rg is calculated from ro by applying a matrix L, describing a right-handed rotation through an 
angle ip about the Z axis of date, rg = Lro- This lag matrix L is 

cos ip sin ^ \ 

-sinV' cosV' . (3.36) 
1/ 

A positive value of ip implies that the peak response on an Earth meridian occurs at a time 6t = ip/u!-E after that 
meridian plane containing ro crosses the tide-producing object, where we is the angular rotation rate of the Earth 
(Eq. 3.163). No significant departures from a zero phase shift have yet been detected: the peak response occurs when 
the meridian plane containing also includes Hp. Some theoretical models predict out-of-phase components due to 
anclasticity. 

The tidal potential at r^ due to the perturbing body at Rp may be expressed as a multipole expansion 
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^) p,{cose) + [^) Ps{cose) 



= U2 + Us, (3.37) 

where only the quadrupole and octupole terms have been retained. Here, G is the gravitational constant, nip is the 
mass of the perturbing body, and P3 are Legendre polynomials, and 9 is the angle between and Rp. While 
the quadrupole displacements are on the order of 500 mm, the mass and distance ratios of the Earth, Moon, and 
Sun limit the octupole terms to a few mm. An estimate of the retardation correction (employing the position of the 
tide-producing mass at a time earlier than that of the observation by an amount equal to the light-travel time) shows 
that this correction is well below 1 mm, and can therefore be neglected. We calculate tidal effects in the terrestrial 
CIO 1903 frame. While is given in this frame, Rp is typically taken from the planetary cphemeris which is defined 
in the SSB frame. The Lorentz transform from the SSB to the geocentric celestial frame given the Earth's orbital 
velocity of 10~''c may be neglected because the maximum solid tide is 500 mm (xlO~^ = 0.05 mm). However, the 
rotation between the celestial frame and the terrestrial frame is large and must be accounted for: 

R-Pt = Q ^R-Po > (3.38) 
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where Q is the rotational transformation between frames defined by Eq. (3.76). The subscripts t and c on Rp indicate 
when the vector is expressed in the terrestrial or celestial frame, respectively. The perturber unit vector Rp has the 
form of the signal propagation vector from the radio source (Eq. 3.3), with the right ascension and declination of the 
perturbing body {up, 5p) replacing the source coordinates (a, 6) and the vector direction reversed: 



— Q ^ 



cos 6p sin ap 
cos Sp cos ap 
sin(5„ 



(3.39) 



The station unit vector takes a similar form in terms of the longitude and geocentric latitude: 



cos (f>s sin Xg 
cos (j)s cos As 
sin 6,^ 



(3.40) 



After the SSB coordinates of the perturbing body are rotated into the terrestrial frame of the station, the scalar^product 
of the unit vectors yields the angle needed to calculate the Legendre polynomials in Eq. (3.37), cos9 = • Rp. 

Once the tidal potential is obtained, the next step is to determine the displacement experienced by a station located 
at Tg. In a local geodetic VEN coordinate system on an elliptical Earth, the solid tidal displacement vector Ssoi is 



(3.41) 



where the Sj'^ {i = 2, 3) are the quadrupole and octupole displacements. The components of Sgoi are obtained from 
the tidal potential as 



(5^ ^ = h,U,/g 



li cos ( 

'dU, 



(3.42) 
(3.43) 

(3.44) 



where hi{i = 2,3) are the radial (quadrupole and octupole) Love numbers, = 2,3) the corresponding tangential 
Love numbers, and As and 4>s are the station longitude and geocentric latitude, and g the surface acceleration due to 
gravity. 



(3.45) 



Values of the various Love numbers range from to 1, since the Earth is neither perfectly rigid nor perfectly elastic. 
For the highest accuracy, distinction must be made between displacements in the radial vs. vertical, and tangential 
vs. horizontal directions, which differ to the extent of the Earth's flattening (about 1 part in 300). 

Reverting to Cartesian coordinates, some algebra then produces the following expressions for the quadrupole and 
octupole components of 6so\ in terms of the coordinates of the station {xs,ys,Zs) and the tide-producing bodies 
(Xp, Yp, Zp): 
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(3.47) 
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(3.48) 
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5(r« • ILpY - 3r^,Rf 



(3.49) 



<5r = a3/5)E 



2RI 



5(r, • R^)^ - rlR^ 



{xsYp - ysXp)/^Jx1+y'l 



(3.50) 
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5(r, ■ Rp)2 - rlRl 



{xsXp + y^Yp) 



where /Up is the ratio of the mass of the disturbing body, p, to the mass of the Earth, and 

R,p — \-^p 1 Yp 1 Zp^ 



(3.51) 



(3.52) 



is the vector from the center of the Earth to that body. The summations arc over tide-producing bodies, of which only 
the Sun and the Moon are normally included. Recent work considers low-level planetary contributions (Hartmann 
and Soffel, 1994; Hartmann and Wenzel, 1994; Williams, 1995), of which Venus' appear to be the most important. 
Mathews et al. (1995) have recently re-examined the basic definitions underlying the derivation of the tidal potential. 
They conclude that the most reasonable definition is one that uses the reference ellipsoid, and thus implies that the 
Love numbers have a slight latitude dependence. 

The above formulation implicitly assumes that the Love numbers hi and li are independent of the frequency of the 
tide-generating potential. A more accurate treatment entails a harmonic expansion of Eqs. (3.46)-(3.51) and use of a 
different set of /ij, U for each frequency component. Each harmonic term is denoted by its historical (Darwin) name, 
if one exists, and the Doodson code (lERS, 1992) (e.g., Ki has Doodson number = 165555). The Doodson notation 
classifies the tidal components according to increasing speed. The correction to the Love number which scales the 
solid tidal radial displacement for the fcth harmonic term at station s is given by 



5hf 



3 ■y/5/247r 5h\ Hk sin (j)s cos (j)s sin(As + 9k) 



(3.53) 



where J/ij is the difference between the nominal quadrupole Love number /12 = 0.609 and the frequency dependent 
Love number (Wahr, 1979), Hk is the amplitude of the kth harmonic term in the tide generating expansion from 
Cartwright and Edden (1973), (ps is the geocentric latitude of the station. A., is the East longitude of the station and 
9k is the kth harmonic tide argument. The Love numbers and tidal amplitudes of the six dominant nearly diurnal 
tides are listed in Table IIL They yield purely vertical station displacements (Naudet, 1994) (in mm) 

(3.54) 
(3.55) 
(3.56) 
(3.57) 
(3.58) 
(3.59) 

where (f>s,^s, and ac are the station geocentric latitude and longitude and Greenwich RA, respectively. The astro- 
nomical arguments F, D, n (mean anomaly of the Sun, mean argument of the latitude of the Moon, mean elongation 
of the Moon from the Sun, and the mean longitude of the ascending lunar node) are defined in Sec. III.D. These dis- 
placements are then summed and the total is used as the first order correction to each station's vertical displacement. 
Horizontal corrections are presently ignored. Note that the largest correction, the Ki term, is identical to that already 
recommended in 1983 by the MERIT standards (Melbourne et al., 1983, 1985). A fairly complete recent treatment 
of the frequency dependence of tidal response has been given by Mathews et al. (1997). With such amendments, the 
solid tide model presented here should be adequate at the mm level. 

To convert the locally referenced displacement, Sso\, which is expressed in the local VEN coordinate system, to 
the Earth-fixed frame, two rotations must be performed. The first, W, rotates by an angle </'s(gd) (station geodetic 
latitude) about the y axis to an equatorial system. The second, V, rotates about the resultant z axis by an angle, 
—As (station longitude), to bring the displacements into the standard geocentric coordinate system. The result is 
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Asol = VW5«oi , 



(3.60) 



where 



' cos ?!>s(gd) 

W=( 1 O '], V=( sinA, cosA, | . (3.61) 

, sin (/)^(gd) 





We calculate geodetic latitudes according to Bowring's (1976) formulation which properly accounts for stations that 
are at some altitude above the reference ellipsoid: 



tan 



el Re sin^e (1 - /) 
^ — ef Re cos^^ 



(3.62) 



where rsp^ = {xg^ +ys'^)^^^ is the station radius from the Earth's spin axis, / (« 1/300) is the geoid flattening factor, 
ef = 2/ — and = ei/(l — el) are the first and second eccentricities squared, and 



= tan 



(1-/) 



(3.63) 



Simpler formulas which assume that the station is on the elUpsoid make errors of order 10~^ radians for a mid-latitude 
station at 1 km altitude. The difference between geodetic and geocentric latitude can affect the tide model by as 
much as / times the tidal effect, or w 1 mm. 



b. Pole tide 

One of the significant secondary tidal effects is the displacement of a station by the elastic response of the Earth's 
crust to shifts in the spin axis orientation. The spin axis is known to describe an approximately circular path of 
« 20-m diameter at the north pole with an irregular period of somewhat less than one year (see Fig. 16 in Sec. VLB). 
Depending on where the spin axis pierces the crust at the instant of a VLBI measurement, the "pole tide" displacement 
will vary from time to time. This effect must be included if centimeter accuracy is desired. 

Yoder (1984) and Wahr (1985) derived an expression for the displacement of a point at geodetic latitude (j)s(gd) and 
longitude Ag due to the pole tide: 



a;|i?E 
'pol = 



sin (j)s{gd) cos 0s(gd) {Px cos Ag + Py sin Xs) hr 
-|-cos20s(gd)(pxCosAs -l-py sinAs) I 

+ sm<j)s(^gi){-Px sin Xs + Py cos Xs) I X . (3.64) 

Here ue is the rotation rate of the Earth, Re the equatorial radius of the Earth, g the acceleration due to gravity 
at the Earth's surface, and h and I the customary Love numbers. Displacements of the instantaneous spin axis 
from the long-term average spin axis along the x and y axes are given by and py. Eq. (3.64) shows how these 
map into station displacements along the unit vectors in the radial (r), latitude (0), and longitude (A) directions. 
With the standard values we ~ 7.292 x 10~^ rad/sec, Re ~ 6378 km, and g w 980.7 cm/sec^, the factor uj'^R/g = 
0.00346. Since the maximum values of px and Py are on the order of 10 meters, and h « 0.6, I ~ 0.08, the maximum 
displacement duo to the pole tide is 1 to 2 cm, depending on the location of the station {(f>s{gd)i ^s)- The locally 
referenced displacement <Jpoi is transformed via the matrix VW of Eqs. (3.61) to give the displacement Apoi in the 
standard geocentric coordinate system. 
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c. Ocean loading 

This section is concerned with another of the secondary tidal effects, namely the elastic response of the Earth's 

crust to ocean tides, which moves the observing stations to the extent of a few cm. Such effects are commonly called 
"ocean loading." The current standard model of ocean loading is general enough to accommodate externally derived 
parameters describing the tide phases and amplitudes at a number of frequencies. The locally referenced three- 
dimensional displacement Socn (components Sj) due to ocean loading is related to the frequencies LUi, amplitudes ^y, 
and phases (pij. In a local Cartesian coordinate system (usually with unit vectors in the vertical, East, and North 
directions) at time t, 

N 

Sj = ^ ^ij cos{uiit + Vi- (fij) . (3.65) 

i=l 

The quantities (frequency of tidal constituent i) and Vi (astronomical argument of constituent i) depend only 
on the ephemeris information (positions of the Sun and Moon). The algorithm of Goad (lERS, 1989; lERS, 1996b) 
is usually used to calculate these two quantities. On the other hand the amplitude and Greenwich phase lag 
ipij of component j are determined by the particular model assumed for the deformation of the Earth. The local 
displacement vector is transformed via Eqs. (3.61) to the displacement Aocn in the standard geocentric frame. 

In present models, the local displacements and their phases, and ipij , arc calculated from ocean tidal loading 
models with as many as 11 frequencies. The eleven components are denoted, in standard notation: K2, S2, M2, 
and N2 (all with approximately 12-hour periods), Ki, Pi, Oi, Qi (24 h), Mf (14 day), Mm (monthly), and Sga 
(semiannual). 

Three ocean loading models have been used over the years. They differ in the displacements that are calculated 
and the number of components that are considered, as well as in the numerical values that they yield for and 
(fij. Schcrncck's results (1983; 1991) are the most complete in the sense of considering both vertical and horizontal 
displacements and all eleven tidal components. They have now been adopted for the lERS standards (1996b). Goad's 
model (1983) was adopted in the MERIT and early lERS standards (1989), but only considers vertical displacements. 
Pagiatakis's (1982) model (Pagiatakis, 1990), based on Pagiatakis et al. (1982), considers only six tidal components 
{S2, M2, N2, i^i. Pi, and Oi). 

An extension of the 1991 Scherneck model (Scherneck, 1993) accounts for modulation of the eleven tidal frequencies 
by multiples of N' , which corresponds to the hmar nodal period (18.6 years). On the assumption that these additional 
terms yield ocean loading amplitudes which are in the same ratio to each main loading term as the companion tides 
are to the main tides, the additional station displacements can be written as 

JV 

^j- = X] X] ^*''^'-J' [('^' + 'nki'^N')t + Vi + UkiN' - (fij)] , (3.66) 

i=l k 

where the k summation extends over all integer multiples Uki of the lunar node A''', and r^i is the ratio of the tidal 
amplitude of each companion k to the tidal amplitude of the parent i. Of 26 such components listed by Cartwright 
and Edden (1973), 20 are estimated to be significant in contributing to the largest ocean loading displacements at the 
0.01 mm level. Table IV shows the multiples nki and amplitude ratios rki for these 20 components. 

In pushing the limits of Earth modeling to below 1 cm accuracy in the mid-1990s, ocean loading station displace- 
ments are one aspect of the models that are undergoing close scrutiny. Initial studies indicate that ocean loading 
amplitudes can be derived from VLBI experiments at an approximate accuracy level of 1-2 mm (Sovers, 1994). When 
estimating parameters, great care must be used in order to avoid singularities due to the identity of certain com- 
ponents of station displacements due to nominally different physical effects ("confounding of parameters"). Since 
some components of ocean loading, solid Earth tides, and ocean tidally induced UTPM variations have the same 
frequencies, certain linear combinations of the station displacements that they cause are identical (see Sees. III.C.2.a 
and III.D.l.c). 

3. Non-tidal station motion 

Many other processes take place in the Earth's atmosphere and in its crust that affect the location of an observ- 
ing station on time scales ranging from seconds to years, and distance scales ranging from local to global. Present 
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knowledge of such processes is somewhat sketchy, but both theoretical and experimental research are starting to 
provide useful results. Local processes include the effects of ground water and snow cover redistribution (Chao, 1996), 
and magma chamber activity in volcanically active areas (Webb et al., 1995). Effects that have more widespread 
repercussions arc atmospheric loading and post-glacial relaxation. In contrast to the tidal effects, whose time depen- 
dence is fixed by the precisely known motions of the bodies of the Solar System, the non-tidal processes do not have 
well-known, periodic time dependencies. The two "global" effects, atmospheric loading and post-glacial relaxation, 
are discussed briefly in the next two sections. Two others have not yet been investigated in detail: non-tidal ocean 
surface height variations (Caspar and Ponte, 1997), and surface displacements caused by pressure variations at the 
core-mantle boundary (Fang et al., 1996). 



a. Atmosphere loading 

By analogy with the ocean tides that were considered in the previous section, a time-varying atmospheric pressure 
distribution can also induce crustal deformation. Rabbel and Schuh (1986) first estimated the effects of atmospheric 
loading on VLBI baseline determinations, and concluded that they may amount to many millimeters of seasonal 
variation. In contrast to ocean tidal effects, analysis of atmospheric loading does not benefit from the presence of a 
well-understood periodic driving force. Otherwise, estimation of atmospheric loading via Green's function techniques 
is analogous to methods used to calculate ocean loading effects. Rabbel and Schuh recommend a simplified form of 
the dependence of the vertical crust displacement on pressure distribution. It involves only the instantaneous pressure 
at the site in question, and an average pressure anomaly over a circular region C of radius R = 2000 km surrounding 
the site. The expression for the vertical displacement (mm) is 

Ar = -0.35Ap - 0.55p , (3.67) 

where Ap is the local pressure anomaly (p — p^,,.^ , mbar) relative to the standard pressure 

Ps^i, = 1013.25 exp(-/i/8.567) . (3.68) 

The pressure anomaly p (mbar) is averaged over the 2000-km circular region mentioned above. The site altitude h 
(km) should be calculated relative to the standard reference ellipsoid, using quantities defined in the solid Earth tide 
section above (Eq. 3.62): 

^"^•^sd ^1 _ el sinVgd 

Note that the reference point for this displacement is the site location at its standard pressure. The locally referenced 
Ar is transformed to the standard geocentric coordinate system via the transformation of Eqs. (3.61). 

Such a rudimentary model has been recently incorporated into VLBI analyses. A mechanism for characterizing the 
pressure anomaly p expresses the two-dimensional surface pressure distribution (relative to the standard pressure) 
surrounding a site as a quadratic polynomial 

p{x, y) = Ap + Aix + A2y + A^x^ + A^xy + Asy^ , (3.70) 

where x and y are the local East and North distances of the point in question from the VLBI site. The pressure 
anomaly may then be evaluated by the simple integration 

P = JJ^dxdyp{x,y) j JJ^dxdy, (3.71) 

giving 

p = Ap+{A3 + A5)Ry4 (3.72) 

and 

Ar = -0.90Ap - 0.14(^3 + A)^^ • (3.73) 
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A quadratic fit to any available area pressure measurements can then determine the coefRcients ^1-5. Future advances 
in understanding the atmosphere-crust elastic interaction can probably be accommodated by adjusting the coefficients 
in Eq. (3.67). As an initial step along these lines, a station-dependent factor / is introduced to scale the local coefficient 
in Eq. (3.73): 

Ar = -0.90(1 + f)Ap - 0.14(^3 + ^5)^^ . (3.74) 

This may account for differing geographical features surrounding different sites. In particular, / may depend on the 
fraction of ocean within the 2000 km radius. Some recent analyses have produced empirical estimates of atmosphere 
loading coefficients for a number of sites (Manabe et al., 1991; van Dam and Herring, 1994; MacMillan and Gipson, 
1994; van Dam et al., 1994). It is not yet clear whether the site-to-site variation of these coefficients is free from other 
systematic errors. 

b. Post-glacial relaxation 

Thick glacial ice sheets covering Scandinavia, Greenland, and Canada melted «10,000 years ago. The removal of 
their considerable weight pressing on the Earth's crust is believed to result in relaxation ("reboimd") that continues 
at present (Tushingham and Peltier, 1991). The magnitude of this motion is estimated to be as large as several 
mm/yr, predominantly in the vertical direction, at sites in and near the locations of ancient glaciers and ice sheets. 
Current theory of deglaciation effects is not yet sufficiently developed to produce definitive results. The parameters 
describing deglaciation history, as well as the rheological properties of the Earth's mantle, are not accurately known. 
Nevertheless, it appears that reasonable parameter choices yield some agreement with empirical measurements of 
baselines in the affected areas (Mitrovica et ai, 1993; Peltier, 1995). 

In summary, models have been presented that describe the four major time-dependent station motions (solid, pole, 
and ocean tides, and atmospheric loading). Each of the locally referenced displacement vectors is then transformed 
to the standard geocentric coordinate system via rotations like Eq. (3.60). After this transformation, the final station 
location in the terrestrial frame is 

rt = ro + Asoi + Apoi -|- Aoc„ -l- Aatm • (3.75) 



D. Earth orientation 

The Earth is approximately an oblate spheroid, spinning in the presence of two massive moving objects, the Sun 
and the Moon. These are positioned such that their time-varying gravitational effects not only produce tides on 
the Earth, but also subject it to torques. In addition, the Earth is covered by a complicated fluid layer, and is not 
perfectly rigid internally. As a result, the orientation of the Earth is a very complicated function of time, which can 
be represented to first order as the composite of a time-varying rotation rate, a wobble, a nutation, and a precession. 
The exchange of angular momentum between the solid Earth and the fluids on its surface, as well as between its crust 
and deeper layers, is not readily predictable, and thus must be continually experimentally determined (Le Mouel 
et al, 1993). Nutation and precession are fairly well modeled theoretically. At the accuracy with which VLBI can 
determine baseline vectors, however, even these models are not completely adequate. 

Currently, the rotational transformation Q of coordinates from the terrestrial frame to the celestial geocentric frame 
is composed of 6 separate transformations (actually 12 rotations, since the nutation, precession, and "perturbation" 
transformations N, P, and ft each consist of 3 rotations) applied to a vector in the terrestrial system: 

Q = J2PNUXY . (3.76) 

In order of appearance in Eq. (3.76), the transformations are the perturbation rotation, precession, nutation, UTl, 
and the x and y components of polar motion. All are discussed in detail in the following four sections. With this 
definition of Q, if rt is a station location expressed in the terrestrial system, e.g., the result of Eq. (3.75), then that 
location expressed in the geocentric celestial system is 

rc = Qr* . (3.77) 
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This particular formulation evolved from the historical development of astrometry, and is couched in that language. 
With modern measurement techniques, such a formulation is esthetically unsatisfactory, but it currently offers a 
practical way for merging VLBI results into the long historical record of astrometric data. The precession and 
nutation matrices (P and N) depend on two angles (nutation in celestial longitude A?/' and in obliquity Ae) that vary 
slowly in the celestial reference frame. On the other hand, the UTPM matrices (U, X, and Y) depend on three angles 
that vary slowly in a terrestrial reference frame co-rotating with the Earth. Thus the desire to have expressions that 
vary slowly with time led to a description that uses 5 angles to specify the orientation of the Earth, rather than the 
absolute minimum of 3 angles. Recently improving measurement accuracy demands high-frequency terms in UXY, 
some of which duplicate parts of N, and thus negate the advantages of the historical decomposition. Esthetically it 
would be much more pleasing to separate Q into two rotation matrices: 



where Q2 are those rotations to which the Earth would be subjected if all external torques were removed (approxi- 
mately UXY above), and where Qi are those rotations arising from external torques (approximately f2PN above). 
Even then, the tidal response of the Earth prevents such a separation from being perfectly realized. Eventually, 
the entire problem of obtaining the matrix Q and the tidal effects on station locations may be solved numerically. 
Note that the matrices appearing in the transformation of Eq. (3.76) are not the same as those historically used in 
astrometry. Since we rotate the Earth from "of date" to J2000.0 (rather than the celestial sphere from J2000.0 to 
coordinates of date), the matrices n, P, and N are transposes of the conventional transformations. We will now 
consider each of the rotations Y, X, U, N, P, and n in detail. 

1. UTl and polar motion 

The Earth's instantaneous spin axis traces a quasi-circular, quasi-periodic path approximately 20 m in diameter 
with a "period" somewhat less than one year, which is known as polar motion (sec Fig. 16 in Sec. VLB). The first 
transformation, Y, is a right-handed rotation about the x axis of the terrestrial frame by an angle 02. Currently, the 
terrestrial frame is the CIO 1903 frame, except that the positive y axis is at 90 degrees east (toward Bangladesh). 
The X axis is coincident with the 1903.0 meridian of Greenwich, and the z axis is the 1903.0 standard pole. 



where ©2 is the y pole position published by lERS. 

The next rotation in sequence is the right-handed rotation (through an angle ©i about the y axis) obtained after 
the previous rotation has been applied: 



In this rotation, ©i is the lERS x pole position. Note that we have incorporated in the matrix definitions the 
transformation from the left-handed system used by lERS to the right-handed system we use. Note also that instead 
of lERS data used as a pole definition, we could instead use any other source of polar motion data provided it was 
represented in a left-handed system. The only effect would be a change in the definition of the terrestrial reference 
system. 

The application of XY to a vector in the terrestrial system of coordinates expresses that vector as it would be 
observed in a coordinate frame whose z axis was along the Earth's ephemeris pole. The third rotation, U, is about 
the resultant z axis obtained by applying XY. It is a rotation through the angle —H, where H is the hour angle 
of the true equinox of date {i.e., the dihedral angle measured westward between the xz plane defined above and the 
meridian plane containing the true equinox of date) : 



Q = Q1Q2 , 



(3.78) 




(3.79) 




(3.80) 




cos H — sin if 
sin H cos H 
1 




(3.81) 
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The equinox of date is the point defined on the celestial equator by the intersection of the mean ecliptic with that 
equator. It is that intersection where the mean ecliptic rises from below the equator to above it (ascending node). 
This angle H is composed of two parts: 

H = hr + a-E , (3.82) 

where hx is the hour angle of the mean equinox of date, and (equation of the equinoxes) is the difference in hour 
angle of the true equinox of date and the mean equinox of date, a difference which is due to the nutation of the Earth. 
This set of definitions is cumbersome and couples nutation effects with Earth rotation. However, in order to provide 
a direct estimate of conventional UTl (universal time) it is necessary to endure this historical approach. 

UTl is defined to be such that the hour angle of the mean equinox of date is given by the following expression 
(Aoki et al, 1982; Kaplan, 1981): 

hr = UTl + 6^ 41°" 50^54841 + 8640184^.812866 T 

+ 0^.093104 - 6'.2 x 10"*^ , (3.83) 

where the quantity 

T = (Julian UTl date - 2451545.0)/36525.0 (3.84) 

represents the number of Julian centuries since J2000.0. An equivalent expression which is normally used is 

hr = 86400"(UT1 Juhan day fraction) + 67310". 54841 

+ 8640184".812866 T + 0^093104 T'^ - 6^2 x 10"*^ . (3.85) 

This expression produces a time UTl which tracks the Greenwich hour angle of the real Sun to within 16™. However, 
it really is sidereal time, modified to fit our intuitive desire to have the Sun directly overhead at noon on the Greenwich 
meridian. Historically, differences of UTl from a uniform measure of time, such as atomic time (UTC or TAI), have 
been used in specifying the orientation of the Earth (see the plot of UTl— TAI in Fig. 17 of Sec. VLB). 

By the very definition of "mean of date" and "true of date" , nutation causes a difference in the hour angles of the 
mean equinox of date and the true equinox of date. This difference, called the "equation of the equinoxes" , is denoted 
by ue and is obtained as follows: 

Q!E = tan-i l^\= tan-i ( ] = ta.n~^ I TT^ h (3-86) 





where the vector 





AT-H (3.87) 



is the unit vector, in true equatorial coordinates of date, toward the mean equinox of date. In mean equatorial 
coordinates of date, this same unit vector is just (1, 0, 0)-^. The matrix N^^ is the inverse (or equally, the transpose) 
of the transformation matrix N, which will be defined below in Eq. (3.101), to effect the transformation from true 
equatorial coordinates of date to mean equatorial coordinates of date. To a very good approximation, 

aE = A'^ cos(e + As) . (3.88) 

The lERS Conventions (1992) had recommended a slightly different expression, which uses the mean obliquity s 
rather than the true obliquity e + As, and can change aE by a few nanoradians: 

aE = A'^cose. (3.89) 

More recently, the lERS (1996b) recommends including two additional terms so that the equation of the equinoxes 
involves the longitude of the lunar node fl (to be defined in Sec. III.D.2): 

aE = AV'cose + 0."00264sinn + 0."000063sin2O . (3.90) 
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This expression should be used in analyses done after 1 January, 1997 (including re-analyscs of old data). That 
particular date is chosen to minimize the discontinuity of q;e (f^ = on 26 February, 1997). One last note on the 
equation of the equinoxes: augmentation of the a priori precession constant will shift the position of the celestial 
ephemeris pole and thus require an adjustment to the equation of the equinoxes 

AaE = ApLscose - AppL , (3.91) 

where ApLS and AppL are the a priori shifts in lunisolar and planetary precession constants, respectively. 

It is convenient to apply UXY as a group. To parts in 10"'^^, XY — YX. However, to the same acciiracy 
UXY ^ XYU. Neglecting terms of 0(8^) (which produce station location errors of approximately 0.006 mm): 

(cosH —sinH — sinOi cosiJ — sin02 siniJ \ 
siniJ cosH - sin 91 sin i? + sin 62 cos if . (3.92) 
sinOi —sin 62 1 / 

Over relatively short time spans. Earth rotation might be modeled as a time-linear function, by analogy with 
tectonic motion of the stations over longer periods. If the x, y components of polar motion and UTl are symbolized 
by Oi-3, and the reference time is to, then this model is 

Q^ = e° + Qi{t-to) , (3.93) 

where 9° are the values of UTPM at the reference epoch. 



a. Tidal UTPM variations 

Tidally induced shifts of mass in the solid Earth, oceans, and atmosphere carry angular momenta which must be 
redistributed in a manner that conserves the total angular momentum. This leads to variations in the orientation and 
rotation rate of the Earth: modification of polar motion and UTl. Such small effects emerged above the detection 
threshold in space geodesy in the early 1990s. Modeling them is important if sub-centimeter accuracy is to be attained 
in the interpretation of VLBI measurements. 

Just as various tidal forces affect the station locations (Sees. III. C. 2. a - III.C.2.c), they also affect polar motion and 
UTl (9i_3). Equations similar to (3.35) may be written for each of the three components of Earth orientation: 

9^ = 9^0 + A9i,oi + A9^oe„ + A9^atm , (3.94) 

where 9j (i=l,3) symbolizes each of the three components of UTPM, 9jo is its value in the absence of tidal effects, 
and the three A terms are the respective contributions of solid Earth, ocean, and atmospheric tides. The next two 
sections describe the current models of solid and ocean tide contributions that are presently included in VLBI analyses. 
At present not enough is known about atmospheric tidal effects, and only loose upper limits can be placed on their 
contribution. 



b. Solid Earth tide UTPM variations 

The pioneering work in tidal effects on Earth orientation was that of (Yoder et al., 1981). It was limited to UTl, 
and included solid tidal and some ocean tidal effects. Their calculated AUTl can be represented as 



N 



AUTl = 



i=l 



Ai sinf kijaj 



(3.95) 



where N=41 is chosen to include all terms with periods from 5 to 35 days. There are no longer-period contributions 
until a period of 90 days is reached. However, these long-period terms {e.g. lERS, 1996b) are already included in the 
results reported by the current Earth-orientation measurenicmt services. The values for fc^ and Ai, along with the 
periods involved, are given in Table 4.37.1 of the Explanatory Supplement (Archinal, 1992a). The aj are just the five 
fundamental arguments defined in Eqs. (3.104 - 3.108) as I, I', F, D, and ft, respectively. 
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c. Ocean tide UTPM variations 



Redistribution of the angular momentum produced by ocean tides affects the Earth's rotation pole position and 
velocity. This effect was first quantified by Brosche (1982), Baader et al. (1983), and Broschc et al. (1989; 1991). 
The dominant effects on polar motion and UTl occur at diurnal, semidiurnal, fortnightly, monthly, and semiannual 
tidal periods. Assuming that the periods longer than one day are adequately accounted for either in modeling which 
combines solid and ocean tidal effects (not strictly true with the Yoder model), or are already present in the a priori 
UTPM series, only the diurnal and semidiurnal frequencies need to be modeled. Further limiting the model to tidal 
components with apparent amplitudes larger than 1 /xs gives eight components. 

For unified notation, again define (Z=l,3) = x, y polar motion and UTl, respectively. Then the ocean tidal 
effects A 6 can be written as 



i=l 



Ail cosi ^ kijaj + UiQiT + tt) j + Bu sin I ^ kijUj + ni{hr + tt) 



(3.96) 



Ail and Bii are the cosine and sine amplitudes that may be calculated from theoretical tidal models (as in the 
work of Brosche) or empirically determined from data (Sovers et al., 1993; Herring and Dong, 1994; Watkins et 
al., 1994; Gipson, 1996). Table V lists the eight dominant terms in the model. These numerical coefficients are taken 
from the empirical results of Sovers et al. (1993), and are representative of a number of recent empirical models. Note 
that the Ki terms contain large arbitrary conventional components, and are thus not solely due to ocean effects. 
Theoretical calculations of polar motion ocean effects have been made by Chao et al. (1996), yielding the amplitudes 
Ail and Bu using models based on TOPEX/Poseidon altimetry (Fu et al., 1994). The theoretical results agree well 
with empirical determinations, and confirm the dominant role of ocean tides in excitation of short-period UTPM 
variations. 

The ocean tidal UTPM effects are also modulated by the 18.6-year lunar node variation N' . As in the case of ocean 
loading station displacements (Sec. III.C.2.c), the contributions A6;' of the companion tides to AQi can be written 
as 



N 



i=l k 



Au cos( J2 hjO^j + "^iihr + tt) + nkiUJN')t + n^N' ) + 



Bu sin{ kijUj + ni{hr + tt) + nkiLON')t + nuN' 



(3.97) 



where the strengths of the companion tides rki are found in Table IV. 

Polar motion as seen in the reference frame B which rotates with the Earth is identical to nutation in the space- fixed 
celestial frame S. A rotational frequency lub = 1 cycle per sidereal day (cpsd) is identical to lus = 0, while lub — 
corresponds to = — 1 cpsd. Generally, nutations at frequency uis correspond to polar motions with frequencies lob 
= — 1 -I- ujs- The retrograde diurnal parts of the polar motion terms with coefficients Ai (12) and Bi (12) corresponding 
to the tidal components listed in Table V (i = 5 to 8) are thus equivalent to components of the nutation model, and 
due care must be taken when both classes of parameters are estimated. 



2. Nutation 

With the completion of the UTl and polar motion transformations, we are left with a station location vector, 
Tdate = UXYrt. This is the station location relative to true equatorial celestial coordinates of date. The last set 
of transformations arc nutation N, precession P, and the perturbation rotation ft, applied in that order. These 
transformations give the station location in celestial equatorial coordinates; 

rc = J^PNrdate • (3.98) 
The first of these transformations, the matrix N, is a product of three separate rotations (Melbourne et al., 1968): 
1. A(e): true equatorial coordinates of date to ecliptic coordinates of date, 
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1 

A(e) = I cose sine | , (3.99) 
— sin £ cos s 

where e is the true obUquity of the ecliptic of date, 

2. C-^(A^): nutation in longitude from ecliptic coordinates of date to mean ecliptic coordinates of date, 

cos Alp sin Alp \ 
C''{AiP)=\ -sinA^ cosAV- , (3.100) 
1/ 

where Aip is the nutation in ecliptic longitude, and 

3. A^(£): ecliptic coordinates of date to mean equatorial coordinates. 

In ecliptic coordinates of date, the mean equinox is at an angle A^. The angle Ae = £ — e is the nutation in 
obliquity, and £ is the mean obliquity (the dihedral angle between the plane of the ecliptic and the mean plane of 

the equator). "Mean" as used in this section implies that the short-period (T < 18.6 years) effects of nutation have 
been removed. Actually, the separation between nutation and precession is rather arbitrary, but historical. The net 
rotation is 

N = A'^(£)C'^(A^/>)A(£) (3.101) 

cos Aip cos £ sin Aip sin e sin Aip 

■cose sin Ai/i cos e cos e cos A^i' + sine sine cos e sin e cos Ai/i — sine cose 
-sine sin Ai/) sin e cos e cos At/; — cose sine sin e sine cos A^/; + cos e cos e^ 

It should again be pointed out (see Sec. III.D above) that this is the reverse of the customary astronomical nutation. 

The 1980 lAU nutation model (Scidelmann, 1982; Kaplan, 1981) is used to obtain the values of Aip and £ — £. The 
mean obliquity is obtained from Lieske et al. (1977) or from Kaplan (1981): 

£ = 23° 26' 21."448 - 46."8150 T - 5."9 x IO'^'T^ + l."813 x lO'^T^ , (3.102) 

T = (Julian TDB date - 2451545. 0)/36525.0 . (3.103) 

The nutation in longitude Atp and in obliquity A£ = £ — £ can be represented by a series expansion of the sines and 
cosines of linear combinations of five fundamental arguments. The latter are nearly linear in time (Kaplan, 1981; 
Hohenkerk et al, 1992) (in radians): 

1. the mean anomaly of the Moon: 

ai=l = 2.35554839 + 8328.69142288T+ 1.5180 x 10"^T^ + 3.1 x 10"^T^ , (3.104) 

2. the mean anomaly of the Sun: 

a2 = l' = 6.24003594 + 628.30195602T - 2.80 x lO'^T^ - 5.8 x lO'^T^ , (3.105) 

3. the mean argument of latitude of the Moon: 

a3 = F = 1.62790193 + 8433.46615832r - 6.4272 x lO'^T^ + 5.3 x lO'^T^ , (3.106) 

4. the mean elongation of the Moon from the Sun: 

04 = D = 5.19846951 + 7771.37714617T - 3.341 x lO'^T^ + 9.2 x lO'^T^ , (3.107) 
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5. the mean longitude of the lunar ascending node on the ecliptic: 



ag = f2 = 2.18243862 - 33.75704593T+ 3.614 x lO'^'T'^ + 3.9 x lO'^T^ 



(3.108) 



With these fundamental arguments, the 1980 lAU nutation quantities can then be represented by 

5 



N 



{Aoj + AijT) sin ^ kjiaiiT) 



and 



N 



(Boj + BijT) cos ( kjiai{T) 



i=l 



(3.109) 



(3.110) 



where the various values of ai, kji, Aj, and Bj are tabulated in Table 3.222.1 of the Explanatory Supplement. 

Inadequacies of the standard nutation model can be corrected by adding various classes of terms to the nutations Atp 
and Ae in Eqs. (3.109) and (3.110). These include the out-of-phase nutations, the free-core nutations (Yoder, 1983) 
with period lu{ (nominally 430 days), and the "nutation tweaks" dtp and 6e, which are arbitrary constant increments 
of the nutation angles Aip and As. Unlike the usual nutation expressions, the tweaks have no time dependence. 
The out-of-phase nutations Atjj° and Ae°, which are not included in the 1980 lAU nutation series, are identical to 
Eqs. (3.109) and (3.110), with the replacements sin ^ cos: 



and 



N 

Ar^Y. 



N 



{A2j + A^jT) cos ( J2 (T) 



{B2j + BsjT) sin^E kjiai{T) 



Expressions similar to these are adopted for the free-core nutations: 

AV-f = (^00 + AiqT) sin(wf T) + {A^o + A30T) cos(u;fT) 



and 



Ae^ = (Boo + BioT) cos(a;fT) + (S20 + ^30^) sin(a;fT) 



(3.111) 



(3.112) 



(3.113) 



(3.114) 



Since the free-core nutation is retrograde, Uf is negative. The nutation model thus contains a total of 856 parameters: 
Aij (j=0,3; j=l,106) and Bij (j=0,3; j=l,106) plus the free-nutation amplitudes A^ (i=0,3), B^ (j=0,3). The only 
nonzero a priori amplitudes are the Aqj, Aij, Bgj, Bij (j=l,106). 

The nutation tweaks are just constant additive factors to the angles Aip and Ae: 



All) ^ Axj} -\- 5tp and Ae ^ Ae -|- Se 



(3.115) 



Deficiencies in the lAU nutation model became clearly evident in the 1980s (Herring et al., 1986). Several methods 
of correcting them are in current use. The first possibility is to use empirically determined values of Sijj, 5e that 
are available from the lERS (see also Fig. 18 in Sec. VLB). With this choice the nutation angles are determined by 
interpolating among results from other VLBl experiments near the date of interest. An alternative is to estimate 5tp 
and 5e directly from the data. 

Another avenue for improving the a priori nutation model is to select one of the published replacements of the 1980 
lAU series. Early work by Zhu et al. (1989; 1990) refined the 1980 lAU theory of nutation both by re-examining the 
underlying Earth model and by incorporating experimental results. Herring (1991) has extended the work of Zhu et al. 
and used geophysical parameters from Mathews et al. (1991) to generate the ZMOA 1990-2 (Zhu, Mathews, Oceans, 
Anelasticity) nutation series. Kinoshita and Souchay (1990) have re-examined the rigid-Earth nutation theory, and 
attempted to include all terms larger than 0.005 mas, in particular planetary terms not present in any previous 



38 



theories. The 263 hmisolar terms have been corrected for the Earth's non-rigidity (Souchay, 1993). Other additions 
and corrections to the Kinoshita-Souchay model are found in Souchay and Kinoshita (1996), Williams (1994; 1995), 
and Hartmann et al. (1996). 

The Kinoshita-Souchay planetary contributions to A'^ and As are 



N 

E 



and 



.10 s .10 

S^j sin (Y^kjiPi{T)\ + C^j cos (^%A(r) 

.10 . .10 s 

Sej sin 



(3.116) 



(3.117) 



where the astronomical arguments are symbohzed by A; the last four /3i are identical with the ai defined above 



{(37 = D = ai, Ps = F = as, /Jg 
planets (in units of radians): 



I = ai, Pio — Q — a^), while the first five are mean heliocentric longitudes of the 



1. 


Venus : 


(3i 


= w = 


3.176146697 + 1021.3285546 T , 


(3.118) 


2. 


Earth : 


/32 


= l-E = 


1.753470314 + 628.30758492 T , 


(3.119) 


3. 


Mars : 


/33 


= Im = 


6.203480913 + 334.06124315 T , 


(3.120) 


4. 


Jupiter : 


/34 


= h = 


0.599546497 + 52.96909651 T , 


(3.121) 


5. 


Saturn : 


P5 


= h = 


0.874016757 + 21.32990954 T , 


(3.122) 



and the sixth is the accumulated general precession: 

l3e=pA= 0.02438175 T + 5.38691 x 10"^ 



(3.123) 



It should be noted that the paper of Kinoshita and Souchay gives expressions for the lunisolar tidal arguments that 
are slightly at variance with the lAU formulas presented above in Eqs. (3.104 - 3.108). These differences may be of 
significance in high-accuracy modeling studies. 

Since the present standard model of nutation is known to be in error by amounts that are large in comparison 
to present measurement capabilities (see Fig. 18 in Sec. VLB), the International Astronomical Union considers it 
important to formulate and adopt an improved nutation model by the end of the century. A working group is 
presently considering variants of the ZMOA and Kinoshita-Souchay models in this connection. In the meantime, the 
lERS (1996a) maintains empirical time series of the nutation angles A'^, As for use in applications that require higher 
accuracy than the 1980 lAU model. 



3. Precession 



The next transformation in going from the terrestrial frame to the celestial frame is the rotation P. This is the 
precession transformation from mean equatorial coordinates of date to the equatorial coordinates of the reference 
epoch {e.g., J2000.0). As was the case with the nutation matrix of Eq. (3.101), this is a rotation whose sense is 
opposite to that of the conventional astrometric precession. It is a composite of three rotations discussed in detail by 
Melbourne et al. (1968) and Lieske et al. (1977): 



cos Z sin Z 
R{-Z) = ( -sinZ cosZ 
1 



(3.124) 



Q(e) 




(3.125) 
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COS C sin C 

R(-0 = I -sinC cosC 0| . (3.126) 
1 

P = R(-C)Q(e)R(-Z) (3.127) 

cos C COS Q cos Z — sm( sin Z cos ( cos 6 sin Z + sin cos Z cos ^ sin 6 ' 
sin ( cos 9 cos Z — cos C sin Z — sin ^ cos 6 sin Z + cos C cos Z — sin ^ sin G 
— sin 6 cos Z — sin G sin Z cos G 

The auxiliary angles Q, Z depend on precession constants, obliquity, and time as 

C = 0.5mr + 0".30188T2 + 0".017998 (3.128) 
Z = 0.5mT + 1".09468T2 + 0".018203 (3.129) 
G = nT - 0".42665 - 0".041833 , (3.130) 

where the speeds of precession in right ascension and declination are, respectively, 

TO = PLS cos £0 - PPL (3.131) 

n=pLssin£o (3.132) 

and pi^s = the lunisolar precession constant (including the geodesic precession), ppi, = the planetary precession 
constant, eq = the obliquity at J2000.0, and T (Eq. 3.103) is the time in centuries past J2000.0. Nominal values at 
J2GGG.0 arc pLS = 5038".7784/cy, ppl = 10".5526/cy; these yield the expressions given by Lieske et al. (1977) and 
Kaplan (1981): 

C = 2306".2181 T + 0".30188 + 0".017998 (3.133) 
G = 2004".3109 T - 0".42665 - 0".041833 (3.134) 
Z = 2306".2181 T + l".09468 + 0".018203 , (3.135) 

Direct estimates of precession corrections can be obtained from observations. The most recent such results (Chariot 
et al, 1995; Walter and Severs, 1996) indicate that the current lAU nominal value of pls is in error by — 3.0±0.2 
milliarcseconds per year. The precession matrix completes the standard model for the orientation of the Earth. 



4. Perturbation rotation 

For exploring departures from the canonical precession, nutation, and UTPM transformations, the standard model 
for the rotation of the Earth as a whole may be modified by small incremental rotations about the resulting axes. 
Define such a perturbation rotation matrix as 

n = A^AyA, , (3.136) 

where 

/I \ 

(3.137) 




with SQx being a small angle rotation about the x axis, in the sense of carrying y into z; 

( 1 -5Qy\ 

Ay = 1 (3.138) 
VQy 1 / 

with 5@y being a small angle rotation about the y axis, in the sense of carrying z into x; and 
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(3.139) 

with 5Qz being a small angle rotation about the z axis, in the sense of carrying x into y. For angles on the order of 
1 arc second we can neglect terms on the order SQ^Re as they give effects on the order of 0.15 mm. Thus, in that 
approximation 

/ 1 5e, -SQyX 
ft = -5Q^ 1 (56^ . (3.140) 

V SQy -se^ 1 / 

Each of the rotation angles can be expressed as a function of time: 

SOi = SGiiT) = 5&^a + Se,T + MT) , (3.141) 

which is the sum of an offset, a time-linear rate, and some higher order or oscillatory terms. In particular, a non-zero 
value of SQy is equivalent to a change in the precession constant, and SQ^ is equivalent to the time rate of change of 
the obliquity s. Setting 

5Q^ = S&y = 6Qz=0 (3.142) 

gives the effect of applying only the standard rotation matrices. 

Starting in the terrestrial frame with the Earth-fixed vector Tq to which tidal effects A are added, we have shown 
(in Sees. III.C through III.D above) how we obtain the same vector r,. expressed in the geocentric celestial frame: 

Fc = nPNUXY(ro + A) . (3.143) 



E. Earth orbital motion 

We now wish to transform these station locations from a geocentric celestial reference frame moving with the Earth 
to a celestial reference frame which is at rest relative to the center of mass of the Solar System. The reason for this 
apparent complication is that the Solar System barycentric (SSB) frame is particularly useful for expressing the source 
positions, formulating gravitational retardation, and performing trajectory calculations of interplanetary spacecraft. 
These SSB station locations will then be used to calculate the geometric delay (see Sec. III. A). We will then transform 
the resulting time interval back to the frame in which the time delay is actually measured by the interferometer — 
the geocentric celestial frame moving with the Earth. 

Let S' be a celestial geocentric frame moving with vector velocity (3c relative to a frame E at rest relative to the 
Solar System center of mass. Due to the Earth's orbital motion, /3 is on the order 10"'*. Further, let r(t) be the 
position of a point {e.g., station location) in space as a function of time, t, as measured in the S (SSB) frame. In 
the E' (geocentric) frame, there is a corresponding position r'(i') as a function of time, t'. We normally observe and 
model r'(i') as shown in Sees. III.C through III.D. However, in order to calculate the geometric delay in the SSB 
frame S, we will need the transformations of r{t) and r'(t'), as well as of t and t', as we shift frames of reference. 
Measuring positions in units of light travel time, we have from Jackson (1975) the Lorentz transformation: 

r'(i') = r(i) + (7 - l)[r(t) • f3]l3/0^ - 7/3* (3.144) 

t' = ^[t-v{t)-/3], (3.145) 

and the inverse transformation: 

r{t) = r'{t') + (7 - l)[r'{t') ■ (3](3/(3^ + 7/3t' (3.146) 

t = ^[t' + r'{t')-^], (3.147) 

where in this section 
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7=(l-/32) . (3.148) 

Let ti represent the time measured in the SSB frame S, at which a wave front crosses antenna 1 at position 
ri(ti). Let r2(ti) be the position of antenna 2 at this same time, as measured in the SSB frame. Also, let t2 be the 
time measured in this frame at which that same wave front intersects station 2. This occurs at the position r2(t2)- 
Following Sec. IILA, Eq. (3.2) we can calculate the geometric delay t2 — ti. Transforming this time interval back to 
the geocentric S' frame, we obtain 

t2' - t[ = 7(i2 - h) - 7[r2(t2) - ri(ti)] • /3 . (3.149) 

Assume further that the motion of station 2 (with barycentric velocity f32) is rectilinear over this time interval. This 
assumption is not strictly true but, as discussed below, the resulting error is much less than 1 mm in calculated delay. 
Thus 



which gives 



and 



r2 {t2 ) = r2{h)+ 13^ {t2 -h) , (3.150) 
r2(t2)-ri(ti) =r2(ii)-ri(ii)+^2(i2-ti) (3.151) 



i2' - i'l = 7(i2 - ti) - 7[r2(ti)-ri(ti)]-/3 - -ffS^ ■ I3[t2 - h] 

= 7(1 - /32 • /3)(i2 - h) - 7[r2(ti) - ri(ti)] • f3 . (3.152) 

This is the expression for the geometric delay that would be observed in the geocentric S' frame in terms of the SSB 
geometric delay and SSB station positions. 

Since our calculation starts with station locations given in the geocentric frame, it is convenient to obtain an 
expression for r2(ii) — ri(fi) in terms of quantities defined in that same geocentric frame. To obtain such an expression 
consider two events [r'i(f'i), r2(ti)] that are geometrically separate, but simultaneous, in the geocentric frame, and 
occurring at time t[. These two events appear in the SSB frame as 

ri(ii) = r[{t[) + (7 - l)[r[{t[) ■ f3]f3/P^ + iPA (3-153) 

and as 

V2{t2) = v'^it'^) + (7- l)[r^(t'i) ./JJ/S/ZJ^ +7/3f; , (3.154) 

where 

t2-h= 7[r^(ii) - ri(t'i)] • /3 • (3-155) 
With these three equations and the expression 

r2(i2)=r2(ti)+/32[i2-ii] (3.156) 

we may obtain the vector r2(ti): 

+j(3t[ - 7/32 [1-2 (<i) - r[{t[)] ■ f3 . (3.157) 
This is the position of station 2 at the time ti as observed in S. From this we obtain: 

r2(ii)-ri(ii) = r^(i;)-r;(i'i) + (7-l)([r^(t'i)-rU4)]-/3)/9//?' 

-7/32[r'2(t'i)-r'i(t'i)]-/3. (3.158) 

As shown in Sec. IILA, the vectors [r2(<i) — ri(ii)] and /32 are all that is needed to obtain t2 — ti for the case of plane 
waves. For curved wave fronts we will need to know the individual station locations in the barycentric frame as well. 
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These we obtain from Eqs. (3.153) and (3.157) with t\ set equal to zero. Setting t'^ = is justified since the origin of 
time is arbitrary when calculating time differences. 

In implementing these transformations, the relationship for the transformation of velocities is also needed. Taking 
differentials of Eqs. (3.146) and (3.147) we have: 

dr = dr' + (7 - l){dr' ■ l3)fi/0^ + ^^dt' , (3.159) 



dt = j{dt' + dr' ■ 13) . (3.160) 
Dividing to obtain dr/dt we obtain for station 2 in the SSB S frame: 

/92 = [/32 + (7 - l)(-92 • P)f3/P^ + 7/3] / [7(1 + 02 ■ 0)] ■ (3.161) 

AH baryccntric positions and velocities required for the calculations in this section are obtained from planetary 
ephemerides in the J2000.0 frame (Standish, 1982; Standish and Newhall, 1985). 
For station 2 relative to the geocentric origin, we have from Eqs. (3.76) and (3.77): 

/32«^^PN^XYr^a;E , (3.162) 

where the inertial rotation rate of the Earth is 

We = 7.292115 x 10"^ radians per second. (3.163) 

This is not a critical number since it is used only for station velocities, or to extrapolate Earth rotation forward for 
very small fractions of a day {i.e., typically less than 1000 seconds). 

The assumption of rectilinear motion can be shown to result in negligible errors. Using the plane wave front 
approximation of Eq. (3.2), an estimate of the error 6t in the calculated delay due to an error A/32 ™ above value 
of 02 is 

ST = k-[r2{h)-r,{t,)]l I —L\ ^rAfS^. (3.164) 

Further, from Eq. (3.161) above (since 7 « 1 + lO^'^), 

A/32 « ^02 ■ (3-165) 

For the vector 02 in a frame rotating with angular velocity w, the error A02 that accumulates in the time interval r 
due to neglecting the rotation of that frame is 

A02 « /32wr . (3.166) 

Thus for typical Earth-fixed baselines, where r < 0.02 s, neglecting the curvilinear motion of station 2 due to the 
rotation of the Earth causes an error of < 4 x 10~^^ s, or 0.012 mm, in the calculation of r. Similarly, neglecting the 
orbital character of the Earth's motion causes a max;imum error on the order of 0.0024 mm. 

A related transformation from the E to E' frames that is needed for models of antenna axis offsets and atmospheric 
effects is an expression for the "aberrated" source direction unit vector Sq- In the SSB S frame, this vector is just the 
negative of the propagation vector of Eq. (3.3): sq = — k. The Earth's diurnal rotational contribution to aberration 
is two orders of magnitude smaller than that due to the Earth's orbital velocity and can be neglected for axis offset 
and tropospheric models. The source unit vector s in the geocentric celestial frame is then given by the Lorentz 
transformation 

s = a^so + a/3/3 (3.167) 

with the coefficients 

and = as [(1 - l)^^ + 7] • (3-168) 



7(1 + ^ -so) ^^Lv^ u ^2 
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Working in a frame at rest with respect to the center of mass of the Solar System causes rclativistic effects due to 
the motion of the Solar System in the "fixed frame" of the extragalactic radio sources to be absorbed into the mean 
positions of the sources and their proper motions. With observational data extending over a sufficiently long time 
span, this motion in inertia! space should be separable from terrestrial and intra-SoIar-Systcm dynamics. The known 
motions of the Solar System barycenter with respect to four standards are listed in Table VI. Both the magnitudes and 
directions are tabulated, with the latter given in terms of Galactic longitude I and latitude b. In order of increasing 
velocity, they are given with respect to the local standard of rest (LSR), Galactic center (GC), Local Group (LG), and 
cosmic microwave background (CMB). Note that the motion of the SSB relative to the CMB is on the order of 10~^c, 
indicating that relativity should play a significant role in modeling its consequences for VLBI. Since the extragalactic 
radio sources are expected to be between the LG and CMB, the motion of the Solar System relative to the "fixed" 
extragalactic reference frame (EGRF) is also expected to be of similar magnitude. Detection of this motion could 
contribute significantly to understanding of the large-scale structure of the universe. 

We note at this point that the extragalactic radio sources have a median redshift z w 1.2, whereas the more 
distant CMB is characterized by redshift z = 10^ (Melchiorri, 1990). If the extragalactic radio sources are assumed to 
be at rest relative to the CMB, one then has an estimate of the SSB velocity in the extragalactic radio frame. Motion 
of the SSB relative to the EGRF potentially contributes to the observed delays in VLBI experiments in three ways: 
geometric, gravitational, and aberrational effects. 

Geometric effects: The geometric effects of galactic rotation can be easily estimated. In the vicinity of the Sun, the 
period for galactic rotation is approximately 240 million years. Thus our angular velocity about the galactic center 
is ~ 27r/2.4 x 10^ = 2.6 x 10~^ radians/year. For sources within the Galaxy, at distances approximately equal to 
our distance from the galactic center, the apparent positions could change by w 26 nrad/yr. An intercontinental 
baseline (10,000 km) could thus be in error by as much as 26 cm/yr (1 nrad « 1 cm) if measurements were based on 
sources within the Galaxy. Since our distance from the galactic center is 8.5 ± 1 kpc ~ 2.7 x 10^ light years (Binney 
and Tremaine, 1987), and most extragalactic radio sources are believed to be w 10^ light years distant, the potential 
baseline error is scaled by the ratio of these two distances, « 3 x 10~^, and becomes only « 0.01 mm/year. Even 
with the present near-20-year history of VLBI data, the purely geometric systematic error due to galactic rotation is 
negligible, and only exceeds the millimeter level for sources closer than 10 million light years. 

Gravitational delay effects: Gravitational delay effects depend on the relative positions of the radio source, the 
intervening masses, and the receiving station. Thus one must consider how the motion of the SSB changes these 
effects. The central bulge of the Milky Way galaxy contains « 10^^ solar masses and it is about 8.5 kpc (w 2 x 10^ 
AU) from the Solar System. The Galactic gravitational delay is thus estimated to be approximately 40 times larger 
than the solar delay. While the solar effect varies on the time scale of a year, the galactic effect varies on a time scale 
of 240 million years. As with the galactic aberration effects, we choose to absorb the static portion (maximum of 4 
arc seconds at 10° from the GC) into the reported radio source positions. We are left to model the variation of this 
effect with time. Given that the SSB's angular motion about the GC is 5 mas/yr and that scintillations from the 
interstellar medium limit VLBI observations at 2 to 8 GHz frequencies to no closer than 10° from the GC, the largest 
expected change in deflection would be < 0.5 /las per year. 

Aberration effects: The 0.0012c SSB CMB velocity will cause large (« 4 arc min) changes in apparent source 
positions due to aberration. To the extent that the SSB velocity does not change on time scales of decades, however, 
one is free to absorb these aberration effects into reported source positions. With such a convention, the problem is 
reduced to considering changes in the SSB velocity. 

Based on the results of studies in galactic dynamics {e.g. Binney and Tremaine, 1987), one may construct a model 
of the SSB's acceleration both toward the Galactic Center (GC) and toward the Local Standard of Rest (LSR). At 
present we have not modeled the motion of the GC relative to the Local Group (LG) nor have we modeled the motion 
of the Local Group relative to the cosmic microwave background. This is equivalent to assuming that the GC-LG 
and LG-CMB velocities are constant on time scales of decades {i.e., acceleration ^ 2 x 10~^^c/yr). 

SSB-LSR velocity: We consider a circular motion of radius ~ 1 kpc with a period of ~ 200 million years. This gives 
a change in velocity 

^SSB-LSR = 2 X 10-12 c/yr (3.169) 

directed toward the local standard of rest. This acceleration will cause a maximum aberration of ~ 0.5 /uas/yr which 
can be ignored because of its small size. 

LSR-GC velocity: We consider a circular motion of radius 8.5 kpc with a period of 240 million years. This gives a 
change in velocity 

KsK-oc = 2 X 10-" c/yr (3.170) 
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directed toward the Galactic center (a = 17*^45™. 5, S = — 28°56'). This velocity will cause a maximum aberration 
of ~ 5 /xas/yr. This aberration effect is at present not included in a priori VLBI models. It is, however, possible to 
estimate the SSB's velocity change from the data. 



F. Source structure 



By analogy with the time dependence of station coordinates caused by tectonic and tidal motion, non-stationarity 
of the source coordinates must also be considered. In VLBI this complication arises because any structure or motion 
of the radio sources is manifested in slight variations of the observed delays. Likewise, in generalization of the model to 
ranging experiments (interplanetary spacecraft, GPS, SLR, LLR) the "structure" of the "source" includes positions 
of the transmitter relative to the center of mass of the satellite, and its orbital path relative to the Earth {i.e., 
the GPS, Lageos, or Moon ephemerides). For VLBI, such effects can be eliminated to a large extent by judicious 
choice of objects in planning the experiments. With the exception of special purpose experiments (such as those 
of Lestrade et al. (1995) on radio stars), the sources are well outside our Galaxy, ensuring minimal proper motion. 
Numerous astrophysical studies during the past two decades have shown that compact extragalactic radio sources 
exhibit structure on a milliarcsecond scale (e.g. Kellermann and Pauliny-Toth, 1981). Such studies are important for 
developing models of the origin of radio omission of these objects. Many radio source structures are found to be quite 
variable with frequency and time (Zensus and Pearson, 1987; Taylor et al, 1994a; Taylor et al., 1994b; Polatidis et 
al., 1995; Thakkar et al., 1995; Henstock et al., 1994; Fey et al., 1996). Survey maps of 187 radio sources in the Taylor 
et al. (1994a) reference showed that the structures of only 8 sources did not exceed a scale of 1 mas at an observing 
frequency of 5 GHz. If extragalactic sources are to serve as reference points in a reference frame that is stable at a 
level below 1 mas, it is important to correct for the effects of their structures in astrometric VLBI observations. 

Corrections for the effects of source internal structures arc; based on work by Thomas (1980), Ulvostad (1988), 
and Chariot (1989; 1990a). A varying non-point-like distribution of the intensity of a source yields time dependent 
corrections to the group delay and phase delay rate observables, Ats and Afg, that may be written in terms of the 
source's intensity distribution 7(s, w, t) as 

At« = d<j>s/div, Afs = d<j>s/dt , (3.171) 

with 



4)s = arctan(-Zg/.Zc) (3.172) 

and 

Z{S}= jjdn 7(s,a;,t){'^^J(27rB.s/A) . (3.173) 

Here (ps is the correction to the phase of the incoming signal, s is a vector from the adopted reference point to 
a point within the source intensity distribution in the plane of the sky, lj and A are the observing frequency and 
wavelength, B the baseline vector, and the integration is over solid angle Jl. Source intensity distribution maps are 
most conveniently parametrized in terms of one of two models: superpositions of Dirac delta functions or Gaussians. 
At a given frequency, the corresponding intensity distributions are written as 

I{s) = Y.^k5{x-xu,y-yk) (3.174) 



or 



I{s) = - — ^exp( - [{x - Xk)cosek + {y- 2/fc)sin6'fe]2/2afc^ 
^ ZirakOk V 

- [{x - Xk)smek -{y- Vk) cos6'fc]V26fc2 ), (3.175) 

where Sk is the flux of component k, and Sk (with components Xk, yu iu the plane of the sky) is its position relative to 
the reference point. For Gaussian distributions, 9k is the angle between the major axis of component k and the ti axis 
(to be defined below), and {ak,bk) are the standard deviations: full widths at half maximum of the (major, minor) 
axes of component k normalized by 2\/21og2. The quantities ■^{»} entering the structure phase (Eq. 3.172) are 
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k 



(3.176) 



for delta functions, and 

^{e} = E exp[-27r2(4[/| + 6^y,2)]| sin _ ^^^^^ 

for Gaussians. Here 

Uk = u cos Ok +v sin 0^ 

Vfe = — wsin^fe +t;cos^fe , (3.178) 

with u and v being the projections of the baseline vector B on the plane of the sky in the E-W and N-S directions, 
respectively. 

Maps may be specified in terms of an arbitrary number of either Gaussian or delta function components. At most, 

six parameters characterize each component: its polar coordinates and flux, and, for a Gaussian, its major and minor 
axes and the position angle of the major axis. The structural correction for phase is computed via Eqs. (3.172) 
and either (3.176) or (3.177). For the bandwidth synthesis delay observable, the structure correction is the slope 
of a straight line fitted to the individual structure phases calculated for each fr(^qucmcy channel used during the 
observation. For example, for Mark 111 data there are typically 8 channels spanning «8.2 to 8.6 GHz at X band, and 
6 channels spanning f«2.2 to 2.3 GHz at S band. Delay rate structure corrections are calculated by differencing the 
structure phases at the two times (Eq. 2.6) used to form the theoretical rate observable. In the case of dual-band 
(S-X) experiments, a linear combination (Eq. 5.11) of the structure corrections calculated independently for each band 
is applied to the dual- frequency observables. 

The practical quc^sticm to he resolved is whether structural corrections based on source intensity distribution maps 
yield significant and detectable corrections to the observables at the present levels of experimental and modeling 
uncertainty. Maps are becoming available for a considerable number of the sources currently observed by VLBI (Fey 
et ai, 1996; available at http://inaia.usno.navy.mil/rorf/rrfid.html). Some of the extended sources show time 
variability on a scale of months. Since the corrections At,, and Afs are quite sensitive to fine details of the structure, 
in such cases new maps may be required on short time scales. Depending on the relative orientation of the source and 
VLBI baseline, the delay correction can be as large as «1 ns for extended sources, which is equivalent to tens of cm. 
Despite the potentially large corrections and difficulties in producing and applying maps, the prognosis appears to be 
good. Chariot (1990b) fomid that data from a multiple baseline geodynamics experiment are adequate to map source 
structures with high angular resolution. The use of such maps for the source 3C 273 yields structure corrections that 
substantially improve modeling of observable delays (Chariot, 1994). Recent work has investigated the problem of 
time dependent structures by interpolating the positions and strengths of the components ejected from the core of 
3C 273, from maps made many months apart, and using such interpolated maps to correct observations made in the 
time period between the maps (Chariot and Sovers, 1998). Since maps of many sources are becoming available from 
imaging studies using the VLBA (Fey et al., 1996), it can be expected that structure corrections will become fairly 
routine for numerous extended radio sources used in astrometric work. 

Empirical evaluation of the eflfects of unknown source structure on VLBI measurements could be made via the 
time rates of change of the source right ascension a and declination 5, based on a time-linear model of the source 
coordinates 

a = ao + a:{t — to) 

6= do + d{t-to) . (3.179) 

Non-zero estimates of the rate parameters a and 5 could arise either from genuine proper motion or from motion of the 
effective source centroid sampled by VLBI measurements. Unambiguous interpretation of such results is problematic, 
but non-zero rates can be used as crude diagnostics for the presence of structure effects. Apparent source position rates 
have been reported by Ma and Shaffer (1991), Jacobs et al. (1993), and Eubanks et al. (1996). Thus far, statistically 
significant rates are not believed to represent true proper motion, but rather to be the consequence of a change in 
the interference pattern caused by spectral redistribution or ejection of components from the central engine of the 
source. Galactic rotation and gravitational radiation (Pyne et al., 1996; Gwinn et al., 1997) may also contribute to 
apparent motions. Apparent position shifts and proper motions may likewise result from gravitational lensing within 
the Galaxy (Hosokawa et al, 1997). 
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In the absence of maps, using parameter estimation from astromctric VLBI data, source structure may also be 
modeled as a superposition of two 6 functions centered at points Pi{xi,yi) and P2{x2, 2/2) respectively, as in Eq. (3.174) 
above. The parameters describing the two components are: 1) flux ratio K = S2/S1, where Sk is the flux of the fcth 
component, 2) component separation s = |s| = IP1P2I, and 3) position angle 6. The position angle is 9 = 0° when 
P1P2 is in the direction of increasing declination S, and 9 = 90° when P1P2 is in the direction of increasing right 
ascension a. Prom Chariot (1990b), the group delay has the following dependence on the structural parameters: 

_ 27TK (1 - K) R[l - cos(2^P)] 

co{l + K) [K^ + 2Kcos{2TrR) + l] ' ^ ' 

where 

R = B-s/X. (3.181) 

Note that this model is not linear in its three parameters. This may complicate parameter estimation in the absence 
of accurate a priori values. For evaluating partial derivatives of r that are needed for parameter estimation, the 
component separation s and baseline B are most conveniently written in terms of their components in the celestial 

system, as 

s = a s sin 9 + d s cos 9 

B = au X + 5 V X . (3.182) 

Then R becomes 

R = s{u sm9 + V cos9) . (3.183) 



G. Antenna structure 

The development in Sees. III. A to III.E outlines how the time delay model would be calculated for two points 
fixed with respect to the Earth's crust. Just as the sources are not point-like (Sec. III.F), the antenna system 
likewise does not necessarily behave as an Earth-fixed point. Not only are there instrumental delays in the sys- 
tem, but portions of the antenna move relative to the Earth. To the extent that instrumental delays are inde- 
pendent of the antenna orientation, they are indistinguishable to the interferometer from clock offsets (which are 
treated in Sec. IV). If necessary, such instrumental delays can be separated from clock properties by careful cal- 
ibration of each antenna system {e.g. Rogers, 1975). Physical motions of each antenna relative to the Earth's 
surface must be considered, however, since they are part of the geometric model. Fig. 9 shows the Deep Space 
Network's Deep Space Station 43, which is a 70-m diameter antenna located in Tidbinbilla, Australia (near Can- 
berra). It can easily be imagined that establishing and measuring the stability of such a large structure is not a 
trivial matter. For comparison, the majority of contemporary steerable antennas are 20 to 35 meters in diameter. 



PICTURE OF AUSTRALIAN DSN ANTENNA 
OMITTED BECAUSE POSTSCRIPT 
UNACCEPTABLE TO LANL ARCHIVE 



FIG. 9. The 70-m diameter DSN antenna near Tidbinbilla, Australia. 
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The model of the antenna structure is divided into several components. First, we consider the geometric delay 
which arises when the two axes used to steer the antenna do not intersect. This effect is then modified to account for 
the slight distortion by tropospheric refraction of the direction in which the antenna is pointed. Variation in ambient 
temperature causes thermal expansion of the structure, and the force of gravity causes deformations that vary as the 
antenna changes orientation. For most types of antennas there is also a subtle effect caused by the rotation of the 
feed horn - relative to a fixed direction on the celestial sphere - as the antenna tracks the sidereal motion of a source. 
A small correction to the model of the zenith troposphere delay is required for antennas with non- intersecting axes. 
Lastly, we discuss the modeling of 'site vectors' which account for experiment-to-experiment variations in the position 
of a mobile antenna relative to a fixed benchmark. 



FIG. 10. Generalized schematic representation of the geometry of a steerable antenna. Two rotation axes and the "reference 
point" P are indicated. 

Before giving the details of our model, let us describe a general antenna pointing system (shown schematically in 
Fig. 10) which applies to all antennas that are steerable along two coordinates. The unit vector s to the aberrated 
source position is shown. Usually, a symmetry axis AD points parallel to s. The point A on the figure also represents 
the end view of an axis which allows rotation in the plane perpendicular to that axis. This axis is offset by some 
distance L from a second rotation axis BE. For many antennas this offset is zero or a few meters, but it can be as 
large as 15 m (for the 43-m diameter antenna at Green Bank, West Virginia). All points on this second rotation axis 
are fixed relative to the Earth. Consequently, any point along that axis is a candidate for the fiducial point which 
terminates this end of the baseline. The point we actually use is the point P. A plane containing the rotation axis 
A and perpendicular to BE intersects BE at the point P. This is a somewhat arbitrary choice, one of conceptual 
convenience. 

Consider the plane W which is perpendicular to the dish symmetry axis, AD, and contains the antenna rotation axis 
A (perpendicular to the plane of Fig. 10). For plane wave fronts this is an isophase plane (it coincides with the wave 
front). For curved wave fronts, however, it deviates from an isophase surface by « L^/{2R), where R is the distance 
to the source, and L is taken as a typical antenna offset AP. For L « 10 meters, R = i?moon = 60i?E ~ 3.6 x 10^ 
m, and the curvature correction L^/(2i?) « 1.4 x 10~^ mm is totally negligible. The source has to be as close as 
R = 50 km, or 0.01i?E, before this deviation approaches a 1 mm contribution to the delay. Consequently, for all 
anticipated applications of radio intcrferomctry using high-gain radio antennas, the curvature of the wave front may 
be neglected in modeling the influence of antenna orientation on the time delay. Likewise, gravitational effects are 
sufficiently constant over a dimension L to permit the use of a single Cartesian frame over the antenna structure, to 
a very good approximation. Provided the instrumental delay of the antenna system is independent of the antenna 
orientation, the recorded signal is at a constant phase delay, independent of antenna orientation, at any point on the 
W plane. Since this delay is indistinguishable from a clock offset, it will be totally absorbed by that portion of our 
model (see Sec. IV). Any delay along the symmetry axis AD up to the position of the feed (inside the dish, along AD 
in the direction s) will likewise be absorbed into the clock model. 

However, for millimeter accuracy one must consider the Earth's orbital velocity of « 10~* c which causes Lorentz 
effects on the order of 10~'*-L. For a 10-m axis offset this amounts to 1 mm. The Lorentz effects from the Earth's 
orbital velocity are accounted for by aberrating the source direction Sq from the SSB frame into the geocentric celestial 
frame, yielding s according to Eq. (3.167). Since the Earth's diurnal rotational velocity is 100 times smaller than the 
orbital velocity, the diurnal Lorentz effects for a 10 m offset would be negligible (< 0.01mm). 
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The advantage of choosing the W plane, rather than some other plane parallel to it, is that the axis A is contained 
in this plane, and axis A is fixed relative to the BE axis by the antenna structure. If I is the length of a line from 
P perpendicular to the W plane, the wave front will reach the Earth- fixed point P at a time Tax = l/c after passing 
through axis A. If tq is the model delay for a wave front to pass between the reference points P of antennas 1 and 2, 
then the model for the observed delay r should be amended as: 

T = To - {Tax2 - Taxi) = TQ + (^1 - k) / C , (3.184) 

where the subscripts refer to antennas 1 and 2. 

To calculate this "axis offset", we follow the treatment given by Wade (1970). First define a unit vector I along 
BE, in the sense of positive away from the Earth. For antennas with their steering axes in altazimuth, equatorial, or 
X-Y mounts, the direction I points toward the local geodetic zenith, the celestial pole, or the local geodetic horizon, 
respectively. Next we define a vector L from P to A, L = | L | • Without much loss of generality in this antenna 
system, we assume that s, L, and I are coplanar. Then: 

L = ±L -h^^ , (3.185) 
|Ix(sxI)| 

where the plus or minus sign is chosen to give L the direction from P to A. When s and L are parallel or antiparallel, 
if the antenna comes closer to the source as L increases, the plus sign is used. Since 

Ix (s X I) = s - I (I • s) , (3.186) 



Z = s-L 



= ±L\Jl- (s- 1)^ , (3.187) 



where the sign choice above is carried through. 

The vector I is first defined in the terrestrial frame for a station s at geodetic latitude ?!>s(gd) and longitude . We 
will consider four standard mount types (Salzberg, 1967): altazimuth mounts, equatorial mounts, and X-Y mounts 
directed either North-South or East- West. Dropping the subscripts s and gd, for an altazimuth mount I is in the 
direction of the local geodetic zenith: 

cos sin A 

I = I cos^cosA I . (3.188) 
sin0 

For an equatorial moimt (also called HA-Dec), I is in the direction of the celestial pole. The ±1 is appropriate for a 
northern/southern hemisphere station with its polar axis pointing toward the North/South pole: 




1=0. (3.189) 



For an X-Y mount antenna, I is in the direction of the lower of the two axes (toward the local geodetic horizon), 
which can be either North-South or East- West depending on the particular orientation of the X-Y axes. For an X-Y 
mount directed North-South: 

- sin sin A 

sin cos A I , (3.190) 



COS0 



while for an X-Y East- West directed mount: 




(3.191) 

For completeness, we mention three unique antennas that have been used in VLBI experiments, but do not fall 
into any of the standard categories considered above. One of these is an antenna that was extensively used by the 
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IRIS project of the National Oceanic and Atmospheric Administration in experiments during the 1980s. It was 
unique because it was an equatorial mount designed for the latitude of Washington, D.C. {4>w), but was deployed 
at Richmond, Florida until it was destroyed in the hurricane of August 1992. The considerable latitude difference, 
and the axis offset of several meters, make it imperative that the antenna geometry be properly modeled. See Severs 
and Jacobs (1996) for details. Two other unique antennas, in Arecibo, Puerto Rico and Nancay, France, are seldom 
used in astrometric and never in geodetic VLBI work. The Arecibo antenna has hardware features which make it 
equivalent to an altazimuth mount. The Nancay array has been treated by Ortega-Molina (1985). 

Since I is given in the terrestrial frame t and s is in the celestial geocentric frame c, we rotate It into the celestial 
frame using the matrix Q given in Sec. III.D: 

Ic = QI* , (3.192) 

where the subscripts t and c indicate the terrestrial CIO 1903 frame and the geocentric celestial frame, respectively. 
With this done, one may now obtain s- Ic and subsequently the axis offset delay from Eq. (3.187). 

Note that for "nearby" sources parallax must also be included {i.e., geographically separate antennas are not 
pointing in the same direction). If Rq is the position of the source as seen from the center of the Earth, and r is the 
position of a station in the same frame, then the position of the source relative to that station is 

R=Ro-r (3.193) 

and in Eq. (3.187) we make the substitution 



(s.I)' = 



(Ro-r)-l '' 
IRo - r| 



(3.194) 



For an extreme antenna offset of 10 m = 10* mm and |r | = i?E = 6.4 x 10^ km, the parallax contribution exceeds 1 mm 
only if the source is nearer than « 6 x lO^'km, or 40% of the Earth-Sim distance. 

Another amendment to the antenna geometry is actually due to the atmosphere. The antenna tracks the apparent 
position of the source after the ray path has been refracted by an angle AE in the Earth's atmosphere, rather than 
along the vacuum source direction vector s (the aberrated source direction geocentric celestial frame). While this 
atmospheric refraction effect is already implicitly included in the tropospheric delay correction through the so-called 
mapping function (Sec. V.B), it must be explicitly accounted for in the antenna axis offset model by modifying s. 
For an extreme case of an elevation angle of 6° the deflection can be as large as 2 x 10""^ radians. Thus, for an 
antenna with dimension L = 10 meters, the component of the antenna model SI « LAE « 20 mm. The Earth-fixed 
source direction vector s is modified to take atmospheric refraction into account on the basis of the change from the 
vacuum elevation E to an apparent value E + AE. In the notation of Sec. V.B, a single homogeneous spherical layer 
approximation yields the bending correction in terms of the dry and wet zenith troposphere delays .^d,w, the first 

oc 

moment of the wet troposphere refractivity /w, = / dqfv,{q) (Sovers and Jacobs, 1996), the dry troposphere scale 



height A, and the Earth radius Re- 

AE = cos-^[cos(£; + ao)/(l + Xo)] - ao , (3.195) 

where 

Xo = {Zd + ^w/Afw)/A , 
ao =cos-i[(l + a')/(l + (T)] , 

(T = A/ Re , 

a' = (jl + a{a + 2)]^/'^ /sinE-l^sin'^E . (3.196) 

This formula agrees with atmospheric ray-tracing results to within 1% at 6° and «15% at 1° elevation. 

A further geometric effect on the antenna structure also has its origin in the environment: variations of the 
temperature cause vertical displacements of the antenna reference point, by analogy with the model for atmospheric 
loading in Sec. III. C. 3. a. For large antennas, these can amount to several mm for ordinary diurnal and seasonal 
temperature variations of 10—20 K. If VLBI data acquired under diverse weather conditions are to be modeled 
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simultaneously, it may be important to account for the vertical motion of the reference point. A rudimentary model 
of the temperature effect assumes that the vertical displacement Ar of the antenna reference point, a distance A/i 
above the ground, is 

Ar = a{T-To)Ah , (3.197) 

where a is the coefficient of thermal expansion and Tq is the reference temperature. The reference temperature is 
taken to be equal to the long-term average temperature at each station. A linear expansion coefficient of 12 ppm per 
kelvin is appropriate for both of the usual constituents of antenna structures, steel and concrete. The vertical motion 
is thus ~ 0.42 mm/K for a 70-m antenna (assumed to have Ah = 35 m). Refinements of the simple model would 
have to consider details of the antenna structure, and allow for thermal lag relative to the ambient temperature {e.g. 
McGinnis, 1977; Nothnagel et al, 1995), including thermal effects on the axis offsets. 

Gravity loading of the flexible dish structure changes an antenna's focal length. Because the component of the 
gravity load along the antenna's primary axis of symmetry is proportional to sinE {E = elevation angle), the changes 
in focal length also have sinusoidal elevation dependence (Clark and Thomsen, 1988). In some antennas the subreflector 
may be moved ( "autofocused" ) to compensate for such gravity deformations and thereby maintain focus. However, 
this procedure does not maintain a constant signal path length through the antenna optics and thus introduces 
systematic errors in the antenna position derived from the measurements unless such changes in path are modeled. 
For example, the Deep Space Network antennas used in weekly Earth orientation measurements (Steppe et al., 1994) 
are designed to be in focus with no subreflector compensation a,t E = 45°. For these 70-m and 34-m high efBciency 
antennas the delays relative to nominal focus a.t E = 45° have been empirically determined to be, in mm (Jacobs and 
Rius, 1989; 1990): 

rsr(70-m) = 77 sin 54 

Tsr (34-m) = 13.5 sin 9.8 , (3.198) 

where the coefficients are known to approximately 5%. This functional form clearly exhibits the relationship between 
subreflector motion and other model parameters. The elevation-dependent term biases the station vertical coordinate, 
while the constant term is equivalent to the clock offset. Other antennas are expected to require similar corrections 
with different coefficients. 

One physical effect that pertains to both instrumental delays and antenna geometry is the differential feed rotation 
for circularly polarized receivers. This is caused by the changing orientation of the antenna feed relative to a fixed 
direction on the celestial sphere. The phase shift 6 is zero for equatorially mounted antennas. For altazimuth mounts, 

tan ^ = cos sin if/ (sin cos (5 — cos sin (5 cos il) , (3.199) 

with (p = station latitude, H = hour angle, and 6 = declination of the source. For X-Y mounts, two cases are 
distinguished: orientation N— S or E— W. The respective rotation angles are 

tan 9 = —sm(p sin H/ (cos <!) cos 6 + sm(p sin 6 cos H) (N — S) 

tan6' = cosiI/(sin(5siniJ) (E-W). (3.200) 

The effect cancels for group delay data, but can be significant for phase delay and delay rate data (up to 100 fs/s for 
the latter). The effect on phase delay is 

r = (02 - ei)/f , (3.201) 

where / is the observing frequency and 9i the phase rotation at station i. 

Another small correction, which couples atmospheric delay and antenna geometry, accounts for the effect of orien- 
tation of hour angle-declination (HA-Dec) and X-Y antennas on the tropospheric path delay. Antennas with non-zero 
axis offsets, whose second rotation axis (A in Fig. 10) moves vertically with changing orientation, have zenith tropo- 
sphere delays that may vary by 1 to 2 mm over the range of available orientations. Equatorial and X-Y mounts fall 
in this class. At low elevation angles this zenith variation is magnified by the mapping function to 1-2 cm. These 
variations must be modeled in experiments whose accuracies are at the millimeter level {e.g., short-baseline phase 
delay measurements). For the highest accuracy, tropospheric mapping functions that depend on altitude also need 
to account for variation of the altitude of the antenna reference point. Reports by Jacobs (1988; 1991) derive the 
corrections based on considering only the dry troposphere component, and include all terms necessary to achieve an 
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accuracy of a few millimeters at the lowest elevations. The correction to be added to the zenith dry tropospheric 

delay is 



SZd = - Zd{L/A) tjj , (3.202) 

where L is the antenna axis offset, A the dry troposphere scale height 8.6 km), and ip is an angular factor that 
varies with the type of mount. For equatorial mounts, 

V' = cos<?icosiI , (3.203) 

where (f> is the geodetic latitude and H the local hour angle east of the meridian. The Richmond antenna correction 
has this form with the Richmond (j) replaced by (j)w for Washington, and i? by a pseudo-hour angle Hr, where 

Hr = arctan^cos E sm{6 — ^ j [cos sin E — sin cos E cos{6 + e)] ^ • (3.204) 

For North-South oriented X-Y mounts, 

V' = sini;/(l-cos2 6'cos2£;)^/2 ^ (3.205) 
where 6 is the azimuth (E of N), and for East-West oriented X-Y mounts, 

tp = smE/{l-sm'^ecos^ Ey/^ . (3.206) 

Finally, we need to consider experiments involving transportable antennas which are placed at slightly different 
locations each time a site is occupied. The most important part of the antenna "motion" between successive site 
occupations is expressed as an offset ("site vector") between the current antenna location and some Earth- fixed 
benchmark. If it is desired to describe a series of such experiments in terms of a single set of site coordinates (those of 
the benchmark), then this offset vector must be determined as part of each observing session. It is usually expressed 
in local geodetic coordinates (vertical. East, and North). Models of this offset vector normally assume that the local 
geodetic vertical direction for the antenna is parallel to that for the benchmark (flat Earth). This implies that the 
changes derived for the benchmark coordinates are identical to those for the antenna coordinates. The error introduced 
by this assumption in a baseline adjustment is approximately cIAB/Re,, where AB is the baseline adjustment from 
its a priori value, d is the separation of the antenna from the benchmark, and Re is the radius of the Earth. To keep 
this error smaller than 1 mm for baselines that differ from a priori values by Ril meter, it is sufficient for d to be < 
6000 meters. 

More troublesome is that an angular error ^6 in determining the local vertical, when using an antenna whose 
reference point is a distance Ah above the groimd, can cause an error of Ah sinSQ « Ah SQ in measuring the baseline 
to the benchmark (Allen, 1982). Unless this error is already absorbed into the measurement of the offset vector, care 
must be taken in setting up the portable antenna so as to minimize ^8. To keep the baseline error < 1mm for an 
antenna height of 10 meters, 50 is required to be < 20 arcscconds. Often plumb bobs are used to locate the antenna 
position relative to a mark on the ground. This mark is, in turn, surveyed to the benchmark. Even the difference in 
geodetic vertical from the vertical defined by the plumb bob may be as large as 1 arc minute, thus causing a potential 
error of 3 mm for antennas of 10 meter height. Consequently, great care must be taken in these measurements, 
particularly if the site is to be repeatedly occupied by portable antennas. 



IV. INSTRUMENTAL DELAY MODELS 

The frequency standards ( "clocks" ) at each of the two antennas are normally independent of each other. Attempts 
are made to synchronize them before an experiment by conventional synchronization techniques such as GPS, but 
these techniques may be accurate to only « 1 /is in epoch and « 10~^^ in rate. More importantly, clocks often exhibit 
"jumps" and instabilities at a level that would greatly degrade interferometer accuracy if not modeled. To account 
for these clock effects, an additional "delay" Tc is included in the model delay, a delay that models the behavior of a 
station clock as a piecewise quadratic function of time throughout an observing session. Usually, however, only the 
linear portion of this model is needed. For each station this clock model is given by 

Tc = Tel + Tc2{t - to) + Tcsit - tof /2 . (4.1) 
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In addition to the effects of the lack of synchronization of clocks between stations, there are other differential 
instrumental effects which may contribute to the observed delay. In general, it is adequate to model these effects as 
if they were "clock-like" . Note that the instrumental effects on delays measured using the multifrequency bandwidth 
synthesis technique (Rogers, 1970) may be different from the instrumental effects on delays obtained from phase 
measurements at a single frequency. This is because the bandwidth synthesis process obtains group delay from the 
slope of phase versus frequency [rgd = {d4>/dLo)\ across multiple frequency segments spanning the receiver passband. 
Thus, any frequency-independent instrumental contribution to the measured interferometer phase has no effect on 
the group delay determined by the bandwidth synthesis technique. However, if the delay is obtained directly from 
the phase measurement ^ at a given frequency w then the phase delay (rpd = 4'/'^) does have that instrumental 
contribution. 

Because of this difference, it is necessary to augment the "clock" model for phase delay measurements: 

Tcpa =rc + ToA{t - to) + Tc5{t - tof/2 , (4.2) 

where Tc is the clock model for bandwidth synthesis observations and is defined in Eq. (4.1). Since present systems 
measure both bandwidth synthesis group delay and phase delay rate, all of the clock parameters described above 
must be used. However, in a perfectly calibrated interferometer, Tc4 = Tcs = 0. This particular model implementation 
allows simultaneous use of delay rate data derived from phase delay, with group delay data derived by means of the 
bandwidth synthesis technique. 

A refinement of the clock model may be required for dual- frequency (S/X) delays. It originates from the differential 
instrumental delay for S- and X-band data, which may be sizeable. For dual-frequency observables, the clock model 
depends on this differential instrumental delay and on the frequencies ws , wx in the individual bands as 

Tc6 - ^s) ■ (4-3) 

The differential instrumental delay Tcq is normally highly correlated with the usual clock offset Td, but under some 
circumstances may convey additional information. For example, if the frequencies u!s,i^x are not exactly the same 
for all measurements, the Tc6 term will not be perfectly absorbed into the clock offset Td- 

To model the interferometer delay on a given baseline, a difference of station clock terms is formed: 

Tc = re(2) - re(l) . (4.4) 

Specification of a reference clock can be postponed until least-squares parameter adjustment, and is of no concern in 
the model description. 



V. ATMOSPHERIC DELAY MODELS 

During its journey from the radio source to the two Earth-based receivers, the radio wave front must pass through 
intergalactic, interstellar, interplanetary, and terrestrial atmospheric media. Both the neutral and charged components 
of these media (neutral molecules and ionized plasma, respectively) modify the propagation speed. This produces 
an overall delay relative to propagation in vacuum. Differential effects among the divergent paths of the two signals 
traversing different regions of the anisotropic medium distort the wave front, and thus contribute to differential delay 
in arrival times at the two Earth-based receivers. Only in the immediate vicinity of the Earth, however, do the two 
ray paths diverge sufficiently to cause measurable differences in arrival times. This section is concerned with models 
which correct the VLBI observables for such additional delays due to propagation effects. It is divided into two parts 
that consider contributions from charged particles and neutral molecules, respectively. In recognition of the dominant 
influence of the near-Earth environment, they are named "Ionosphere" and "Troposphere" . 



A. Ionosphere 

For a medium composed of charged particles (plasma), Spitzer (1962) gives the refractive index at frequency v in 
the quasi-longitudinal approximation: 



n = 



1 



1/2 



(5.1) 
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where 9 is the angle between the magnetic field B and the direction of propagation of the wave front. The plasma 
frequency Vp is 



/ N 1/2 

fpcVo/7rj Ri8.98p^/^ (Hz) 



and the electron gyrofrequency i^g is 



eB 
27rmc 



2.80 X lO^B (Hz) 



(5.2) 



(5.3) 



Here p is the number density of the electrons (m^''), c the speed of light, ro the classical electron radius, e and m the 
electron charge and mass, and B is measured in gauss. 

Tables VII and VIII give approximate values of the plasma frequency Up and gyrofrequency Vg for the three regimes 
along a radio signal's ray path: Earth, interplanetary, and interstellar space. The plasma parameters p and B 
represent typical conditions in these three regions, with the upper limit adopted in cases of substantial variability 
(Lang, 1980; Zombeck, 1990). For example, the "Earth" values are appropriate to the daytime ionosphere. For ray 
paths close to the Sun, however, the tabulated "interplanetary" values are several orders of magnitude too small. 
Intergalactic electron densities and magnetic fields are at least two orders of magnitude smaller than the tabulated 
interstellar values. Given the S-band {vs = 2.3 GHz) and X-band {fx = 8.4 GHz) frequencies typically used in geodetic 
and astrometric VLBI experiments, the corrections for plasma and gyrofrequency effects can be parametrized by the 
ratios of and Ug to i^s.x respectively. It can be seen from Table VII that there is an order-of-magnitude falloff in 
the plasma frequency as we proceed outward from the Earth to interplanetary and interstellar regions. Table VIII 
shows that electron gyrofrequency effects are a factor of 10 smaller than plasma effects near the Earth, and negligible 
in interplanetary and interstellar environments. 

Relative to a perfect vacuum as a reference, the contribution Ap^ to the VLBI phase delay Tpd for a monochromatic 
signal traversing a medium of refractive index n is 



.,. = l/,„-:,.,.-l/(^)'[i.l(i)%l(^)%.. 



dl 



where 



For 8.4 GHz, we may approximate this effect to parts in 10® — lO'' by: 



1/2 



^Pd ~ -7:2 
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(5.4) 



(5.5) 



(5.6) 



(5.7) 



and where Ig is the total number of electrons per unit area along the integrated line of sight. The angular brackets 
symbolize a geometrical average. If we also neglect the term {vg cos Q) / v, then the expression for Ap^ becomes simple 
and independent of the direction of travel of the wave front through the plasma: 



A. 



pd 



{5. 



This additional delay is negative. Thus, a phase advance occurs for a monochromatic signal. Since phase delay 
is obtained at a single frequency, observables derived from the phase delay {e.g., phase delay rates) experience an 
increment which is negative (the observable with the medium present is smaller than it would be without the medium). 
In contrast, group delays measured by a technique such as bandwidth synthesis [rgd = (9(/)/9cj)] experience an additive 
delay which can be derived from Eq. (5.8) by differentiating <p = vAp^ with respect to frequency: 



Agd = q/iy'^ 



(5.9) 
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The group delay is of the same magnitude as the phase delay advance. For group delay measurements, the measured 
delay is larger with the medium present than without the medium. 

The uppermost component of the Earth's atmosphere, the ionosphere, is a layer of plasma whose density peaks at 
about 350 km altitude, but is widely distributed between 80 and 1000 km above the Earth's surface. It is created 
primarily by the ultraviolet portion of the sunlight and the solar wind, and is a superposition of three Chapman 
layers (Bassiri and Hajj, 1993): E, Fi, and F2, with height parameters of 110, 210, and 350 km, respectively. During 
daytime in a year near the maximum of the 11-year solar activity cycle the electron density peaks in the F2 layer at 
about 3.7 X 10^^ electrons/m"^. For a typical ionosphere, Ag^ 0.1 to 2 ns at local zenith for v = 8.4 GHz. This 
effect has a maximum at approximately 1400 hours local time and a broad minimum during local night. For long 
baselines, the effects at each station can be quite different. Thus, the differential effect may be of the same order as 
the maximum. 

For the interplanetary medium at an observing frequency of 8.4 GHz and assuming that the ionized region extends 

over lO'' km, a single ray path experiences a delay of approximately 6x10 ^ s in traversing the Solar System. However, 
the differential delay between the ray paths to the two stations on the Earth is considerably smaller, since the gradient 
between the two ray paths should also be inversely proportional to the dimensions of the plasma region {i.e., millions 
as opposed to a few thousand kilometers). An interstellar signal from a source at a distance of 1 megaparsec (3 x 10"'^^ 
km) would experience an integrated plasma delay of approximately 600 seconds at a frequency of 8.4 GHz if the 
plasma density were that of Table VII everywhere along the ray path. In the near vacuum of intergalactic space, the 
average plasma density is « 2 orders of magnitude lower than this, reducing the delay to «10 s. In this case, however, 
the typical dimension is also that much greater, and so the differential effect on two ray paths separated by one Earth 
radius is still not as large as the differential delay caused by the Earth's ionosphere. 

Plasma effects can best be calibrated by the technique of observing the sources at two frequencies, vi and V2, where 
1^1,2 ^ Vp and where \u2 — vi\l{v2 + vi) ~ 1- For the VLBI delays t at the two frequencies vi and 1^2 we obtain 

Tvi=T + q/yl and Tv2 = T + q/vl. (5.10) 

Multiplying each expression by the square of the frequency involved and subtracting, 

T = aTv2+b Tv\ , (5-11) 

where a = 1^2/(^2 ~ ^1) ^^'^ ^ — ~^i/(^2 ~ ^i)- This linear combination of the observables at two frequencies thus 
removes the charged particle contribution to the delay. For uncorrelated errors in the frequency windows, the overall 
error in the derived delay can be modeled as 

<7^=aV^,+feV2^^ . (5.12) 

Modeling of other error types is more difficult and will not be treated here. Since the values of a and b are independent 
of q, these same coefficients apply both to group delay and to phase delay. 

If we had not neglected the effect of the electron gyrofrequency in the ionosphere, then instead of Eq. (5.11) above, 
the combined S/X delay would have been 

T = aT„2 + bT„i + q {ug cos6)/[f2Z^i(f2 - I'l)] , (5.13) 

where a and b are defined as in Eq. (5.11). If the third term on the right-hand side is expressed in units of the 
contribution of the ionosphere at frequency f2, we obtain 

r = aT„2 + bTui + Apdf2 {vg cos Q)/[vi{v2 + ^1)] ■ (5-14) 

For X band Apj 0.1 to 2 ns at the zenith. When using S band as the other frequency in the pair, this third term 

is « 2 X 10~^Apd cosQ « 0.02 to 0.4 ps at zenith. In the worst case of high ionospheric electron content, and at low 
elevation angles, this effect could reach 1 mm of total error in determining the total delay using the simple formula of 
Eq. (5.11) above. Notice that the effect becomes much more significant at lower frequencies. Bassiri and Hajj (1993) 
present a more detailed discussion of higher order ionospheric effects. 

To avoid the complications of modeling the ionosphere, most contemporary VLBI experiments are performed at 
two frequencies. The dual-frequency calibration outlined above then permits us to consider the ionosphere problem 
to be solved at currently required accuracy levels, and to ignore ionospheric propagation effects in modeling. For 
completeness, we present a model of ionospheric delay that can be used for single-frequency experiments. In these 
cases, the interferometer model must use whatever measurements of the total electron content are available. The 
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model is very simple: the ionosphere is modeled as a spherical shell whose lower boundary is at the height hi above 
the geodetic surface of the Earth, and the upper boundary is at the height /12 above that same surface (see Fig. 11). 
For each station the ionospheric delay assumed to be 



where 



T, = kghS{E)/y^ 



k = 0.1cro/27r . 



(5.15) 



(5.16) 



/(, is the total electron content at zenith (in electrons per meter squared xlO""'^^), and g = 1(— 1) for group (phase) 
delay. E is the apparent geodetic elevation angle of the source, S{E) is a slant range factor discussed below, and 1^ is 
the observing frequency in gigahertz. 
The slant range factor (see Fig. 11) is 



S{E) 



IR? sin^ E + 2Rh-2 



hi 



i?2 sin^ E + 2Rhi 




hi) 



(5.17) 



which is strictly correct for a spherical Earth of radius R, and a source at apparent elevation angle E. The model uses 

this expression, with the local radius of curvature of the geoid surface at the receiving station R taken to be equal to 
the distance from the station to the center of the Earth. The model also assumes this same value of R can be used 
at the ionospheric penetration points, i.e., 



Ri = R + hi 



(5.18) 




EARTH CENTER 



FIG. 11. Geometry of the spherical ionospheric shell used for ionospheric corrections. The slant range through the shell is 
I2 — h elevation angle E. 

This is not strictly true, but is a very close approximation, especially in view of the relatively crude nature of the 
total electron content determinations on which the model also depends. The total ionospheric contribution on a given 
baseline is 



Tr=T,(2)-T,(l) , 



(5.19) 



where the arguments 1 and 2 identify the stations. The ionospheric total electron content le is obtained by some 
external set of measurements such as Faraday rotation or GPS techniques. Such external measurements, in general, 
are not along directions in the ionosphere coincident with the ray paths to the interferometer. Thus, for each antenna. 
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it is necessary to map a measurement made along one ray path to the ray paths used by the interferometer. Many 
different techniques to do this mapping have been suggested and tried; Wilson et al. (1995) discuss recent progress. 

The deficiencies of these ionosphere models for single-frequency observations are compounded by the lens effect of 
the solar plasma. In effect, the Solar System is a spherical plasma lens which causes the apparent positions of the 
radio sources to be shifted from their actual positions by an amount which depends on the solar weather and on the 
Sun-Earth-source angle. Since both the solar weather and the Sun-Earth-source angle change throughout the year, 
very accurate single-frequency observations over the time scale of a year are virtually impossible unless simultaneous 
auxiliary experiments are performed. 

B. Troposphere 

The lower few tens of kilometers of the Earth's atmosphere are known as the troposphere. In contrast to the 
ionosphere, this layer is to a good approximation electrically neutral (Houghton, 1986). Radio signals passing through 
the troposphere experience delay, bending, and attenuation relative to an equivalent path through a vacuum, because 
the index of refraction is not equal to the vacuum value. The additional delay is « 2 m at zenith and increases to « 20 m 
at 6° above the horizon; bending amounts to « 0.1° at elevations 6° above the horizon. This makes it imperative for 
accurate VLBI models to account for tropospheric delay. We review the refractivity of the moist air composing the 
troposphere, discuss mapping functions which model the integrated path length through the troposphere, and finally 
consider the limitations to the mapping function models due to azimuthal asymmetry and turbulence in the water 
vapor distribution. 

Permanent and induced dipole moments of the molecular species present in the atmosphere modify its index of 
refraction and thus delay the passage of radiation at microwave frequencies. The excess path delay caused by the 
troposphere relative to a vacuum is 



where n is the index of refraction and S represents the signal's path through the troposphere. Since departures of n 
from unity are small, normally the "refractivity" N = 10^(n — 1) is used instead of the index of refraction. Detailed 
discussions of the refractivity of moist air are found in Thompson, Moran and Swenson (1986, Chap. 13), Bean and 
Dutton (1966), or Thayer (1974). To summarize these reviews, the refractivity of moist air at microwave frequencies 
depends on the permanent and induced dipole moments of the molecular species present in the atmosphere. The 
principal species, nitrogen and oxygen, have no permanent dipole moments, and contribute only via their induced 
moments. On the other hand, water vapor docs have a substantial dipole moment. Induced and permanent dipole 
moments contribute to the refractivity as oc p/T and oc p/T^, respectively (Debye, 1929), where p and T are the 
pressure and temperature of the species under consideration. This is the basis of the Smith- Weintraub equation 
(Smith and Weintraub, 1953): 



where p^ and Py are the partial pressures of "dry" air and water vapor in units of millibars. The coefficients are taken 
from Thayer (1974). The first term represents the aggregate induced dipole refractivity from all the dry constituents, 
and the second and third terms the induced dipole and permanent dipole refractivity from water vapor, respectively. 
Thus the troposphere problem becomes a matter of modeling pressure and temperature {p^ , Pv, T) along the ray 
path. Typically this problem is simplified by assuming the troposphere to be a mixture of perfect gases with unit 
compressibilities which are nearly in hydrostatic equilibrium. Parenthctic;ally, as noted by Thompson, Moran and 
Swenson (1986), at visible wavelengths water vapor has a smaller influence on refractivity, so that the third term is 
not needed to model optical refractivity. This simplifies troposphere modeling for laser techniques. 




(5.20) 



N = 77.6 ^ + 64.8 ^ + 3.776 x 10' 



,5 Pv 

rp2 



(5.21) 
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FIG. 12. Mixing ratios of four major constituents of tlie Earth's atmosphere vs. altitude (based on Goody and Yung, 1989). 

One of the biggest obstacles to modeling the tropospheric delay is the non-uniform and highly variable distribution 
of water vapor. In contrast to the smooth exponential decrease of N2 and O2 concentrations with altitude, the H2O 
concentration generally decreases more rapidly within 1-2 km of the surface, and normally becomes negligible above 
~5 km altitude (Arya, 1988; AGU, 1995). Figure 12 shows the altitude dependence of the mixing ratios of four 
major atmospheric constituents. The fractions of O2, Ar, and CO2 are seen to remain nearly perfectly constant from 
the surface up to 50 km, while the water vapor fraction falls by over 3 orders of magnitude up to the tropopause 
altitude (wlS km) before leveling off at 7 ppm. Water vapor content can also easily vary by 50% during the few hours 
around sunrise and sunset. Weather changes have more pronounced effects on the wet than on the dry troposphere. 
The nitrogen and oxygen concentrations essentially scale with pressure while water vapor undergoes much larger 
concentration changes and redistribution with changing weather. 

The next step in generating the troposphere model is to obtain the integrated refractivity for a signal propagating 
from the zenith. The dry portion, primarily oxygen and nitrogen, is very nearly in hydrostatic equilibrium. As a 
result the zenith delay does not depend on the details of and T along the signal path. The delay can be quite 
accurately estimated simply by measuring the barometric pressure at the surface. The dry zenith delay (m) is 
related to the total surface pressure p (mbar) as 

Zd = 2.2768 X 10"V(l - 0.00266 cos2(/) - 0.00028/i) , (5.22) 

where the factor in the denominator corrects for a non-spherical Earth (Saastamoinen, 1972), (j) is the station latitude, 
and h is the station altitude in km (see Eq. 3.69). Davis et al. (1985) discuss the reason for using the total surface 
pressure p instead of the dry partial pressure p-^ . They also explain the related and subtle distinction between 'dry' 
and 'hydrostatic' delay, which has been glossed over in the present discussion. Typically, at sea level in the local 
zenith direction, the additional dry (or hydrostatic) delay that the incoming signal experiences due to the troposphere 
is « 7.7 nanoseconds or 2.3 meters. 

At each station i the delay experienced by the incoming signal due to the troposphere can most simply be modeled 
using a spherical-shell troposphere consisting of wet and dry components Twi and Td,: 

Ttr(i) = Twi + Tdi . (5.23) 

The total troposphere delay model for a given baseline is then: 

Ttr = Ttr(2) - Ttr(l) • (5.24) 

If Ei is the unrefracted geodetic elevation angle of the observed source at station i, we have (dropping the subscript i): 

Ttr{E) = Md{E)Zd + M^{E)Z^ , (5.25) 

where Z<j,w is the additional (dry, wet) delay at local zenith due to the presence of the troposphere, and A4 are the 

so-called "mapping functions" which relate delays at an arbitrary elevation angle E to the zenith delays. Note that 
since E is the unrefracted elevation angle, refraction effects are contained within the mapping functions A^d,w These 
will be considered in some detail below. 
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For some geodetic experiments, the observed delay can be accurately calibrated for the total tropospheric delays at 
the two stations, which are in turn calculated on the basis of surface pressure measurements for the dry component, 
and water- vapor radiometer (WVR) measurements for the wet component. At a fixed location, the dry zenith delay 
(m) is related to the surface pressure p (mbar) via Eq. (5.22). The wet zenith delay Z„ can be inferred from WVR 
measurements performed in the vicinity of the VLBI stations at the time of the experiment. These corrections can 
be removed and replaced by an alternate model if desired. In the absence of such external calibrations, it was found 
that estimating the zenith delay as a linear function of time can improve troposphere modeling considerably. The dry 
and wet zenith parameters are written as 

Zd,w = <w + ZdAi - io) , (5.26) 
where to is a reference time. The time rates of change Zd,w may then be estimated from fits to the data. 



1. Mapping functions 

The term "mapping function" is used to describe the relation between the tropospheric delay at zenith and an arbi- 
trary angle E above the horizon. Throughout the history of VLBI, extensive attention has been paid to tropospheric 
mapping functions, in view of the dominance of tropospheric delay mismodeling in the error budget. The simplest 
way to relate the tropospheric delay for an oblique path to the delay for a signal received from directly overhead is 
to assume a flat Earth covered by an azimuthally symmetric troposphere layer. Since the delay is proportional to 
the path length through this layer, the mapping function is then simply equal to the cosecant of the elevation angle: 
M{E) = 1/ sinE. Even in the early days of VLBI it was found that this approximation was inadequate, and Marini 
(1972) considered corrections accounting for the Earth's curvature. This led to a mapping function in the form of a 
continued fraction, 

M = ^- . (5.27) 

sm E H ; 

n 

shiE ■ 
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The simplest function, which was widely used in early VLBI modeling in the 1970s and early 1980s, was obtained 
by C. C. Chao (1974). It truncates Marini's continued fraction and replaces the second sinE by tani? in order to 
force = 1 at zenith: 

Md,w = . (5.28) 

sm£; + — 

tan E + Bd,w 

The dry and wet coefficients were determined from empirical fits to ray tracing results through measured average 
atmospheres: = 0.00143, = 0.0445, = 0.00035, and = 0.017. It was claimed to be accurate at the level 
of 1% at = 6°, and to become rapidly more accurate as zenith is approached. 

As more accurate measurements in the 1980s demanded more accurate troposphere modeling, numerous improved 
mapping functions were developed. Many of these have mathematical forms that are slight variants of Marini's con- 
tinued fraction, and contain constants derived from analytic fits to ray-tracing results either for standard atmospheres 
or for observed atmospheric profiles based on radiosonde measurements. The functions of Davis et al. (1985), Ifadis 
(1986), Herring (1992), and NieU (1996) (NMF, Niell Mapping Function) faU into this category. Most of them contain 
parameters that are to be determined from surface meteorological measurements. The NMF function is unique in 
that it attempts to represent global weather variations analytically as a function of location and time of year, and 
contains no adjustable parameters. The NMF mapping functions provide an excellent model for situations in which 
data are not available to supply optimum parameters for some of the alternative mappings. Their seasonal, latitude, 
and altitude variations arc based on interpolations of the U. S. standard atmospheres of Cole et al. (1965). Niell 
presents evidence that short term (hours to days) surface temperature variations are unimportant compared to the 
seasonal, latitude, and altitude variations of the temperature in the region from a few km up to the tropopause. 
Another tropospheric mapping function, due to Lanyi (1984) is unique in that it does not fully separate the dry and 
wet components and thus attempts to give a more faithful representation of the physical effects. It also contains the 
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most complete set of atmospheric parameters, and we will therefore present some details for this function. Several re- 
views have recently evaluated the multitude of tropospheric mapping functions that are now available: Gallini (1994), 
Mendes and Langley (1994), and Estefan and Sovers (1994). 

Because of its potential to provide the most complete description of the tropospheric delay, including use of in situ 
temperature profile data, we will give some details of the Lanyi mapping function. Motivation for and details of its 
development were given by Lanyi (1984). It is based on an ideal model atmosphere whose temperature is constant 
from the surface to the top of the inversion layer hi, then decreases linearly with height at a rate W (lapse rate) from 
hi to the tropopause height /12, and is constant again above /i2- Lanyi's approach was to develop a semi-analytic 
approximation to the ray trace integral (5.20) which retains explicit temperature vs. altitude profile parameters that 
can be supplied from meteorological measurements. The goal was to provide an accurate approximation to ray tracing 
in a form that was less computationally intensive than performing a full numerical ray trace for every observation. 
The mapping function is expanded as a second order polynomial in and (plus the largest third order term). 
Unlike all the other functions mentioned above, the Lanyi mapping function is nonlinear in Z^, Z„. It differs from 
other mapping functions in the sense that it does not conform to the linear expression (5.25). In particular, we note 
that there is a Z^Z^ cross term which couples the dry and wet delays. This nonlinear term is present because for 
a given dry (wet) zenith delay the geometric path through the atmosphere is not independent of the amount of wet 
(dry) zenith delay. Hence the coefficients of the nonlinear terms have a subscript "6" to indicate that these terms 
arise from the bending of the signal path through the atmosphere. 

Here we give only a brief summary of the functional form. The tropospheric delay is written as: 



The quantities Z^ and Z^ are the zenith dry and wet tropospheric delays, while A is the dry atmospheric scale height, 
A = kTo/mgc, k = Boltzmann's constant, Tq = daily average surface temperature, m = mean molecular mass of dry 
air, and gc = gravitational acceleration at the center of gravity of the air column. With the standard values of k, m, 
gc = 978.37 cm/s^, and an average mid-latitude temperature Tq = 292 K, the scale height A = 8.6 km. The dry, wet, 
and bending contributions to the delay, Fd{E), F^{E), and FhiM2M3Mi^)^ ^^'^ expressed in terms of moments of the 
refractivity. The latter are evaluated for the ideal model atmosphere and thus give the dependence of the tropospheric 
delay on the four model parameters Tq, W, hi, and /i2- Note that Lanyi's formulation (Eq. 5.30) differs from the 
simple model (Eq. 5.25) in the presence of the "bending" terms Fbi-4. These account for the influence of the dry and 
wet constituents in bending the incoming ray path. 

Four meteorological parameters describe the temperature vs. altitude profile in the Lanyi model. These have already 
been mentioned above: the surface temperature Tg, temperature lapse rate W, inversion and tropopause altitudes hi 
and Table IX summarizes their standard values, and also gives the approximate; sensitivities of the tropospheric 
delay to the meteorological parameters. These values are calculated at 6° elevation, which is the approximate lower 
limit of validity of the Lanyi model. At this elevation, the ray path traverses the equivalent of approximately 10 "air 
masses". A fifth parameter, the surface pressure po, may be used to calibrate the dry zenith delay via Eq. (5.22). 
The full potential of the Lanyi function can only be realized if complete meteorological information is available for 
the time and place of a VLBI experiment. When interpolated standard global values (e.g. Cole et al., 1965) are used 
for the four parameters, it is essentially equivalent to the NMF mapping. 

While an exceptional amount of research has been devoted to improvement of tropospheric mapping functions 
during the past two decades, it is nevertheless becoming obvious that present measurement accuracy also demands 
characterization of azimuthal asymmetry and short-term weather variations. This may be achievable via real-time 
auxiliary measurements (such as WVRs), or statistical models of temporal and spatial correlations. 

2. Limitations of mapping 

In recent years the limitations of the mapping function have become apparent in VLBI observations. Winds at 
high altitudes, unusually strong lee waves behind mountains {e.g., Owens Valley, California), and very high pressure 
gradients may all limit the accuracy of an azimuthally symmetric dry troposphere model based on measurements 



Ttr = F{E)/ sin E , 



(5.29) 



in order to factor out the l/sin£ "flat Earth" model, where 



F{E) = Fd{E)Zd+F^{E)Z^ 

+ [Fi,i{E)Zd^ + 2Fb2{E)ZAZ^ + Fb3{E)Z^^yA + FM{E)Zd^/A'' . 



(5.30) 
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of surface barometric pressure. Rough estimates indicate that, except in such umisual cases, errors in the simple 
model cause sub-centimeter errors in the baseline. The accuracy in more typical situations is expected to be limited 
by horizontal refractivity gradients caused by equator-to-pole (North-South) gradients in temperature, pressure, and 
humidity. These gradients were first accoimtcd for in analyses of satellite laser ranging experiments by Gardner (1976; 
1977). East-West gradients caused by motion of weather systems passing over a site may also contribute to the 
breakdown of the simplified mapping function model. The limits of validity of the azimuthal symmetry assumption in 
VLBI analysis arc starting to be investigated (MacMillan, 1995; MacMillan and Ma, 1997; Chen and Herring, 1997). 

While the limitations of dry mapping are small and in the case of North-South gradients may be modelable, the 
wet component of the atmosphere (both water vapor and condensed water in the form of clouds) is not so easily 
modeled. It is known to be highly variable in time and space (Flcaglc and Busingcr, 1980). The experimental 
evidence (Resch, 1984) is that it is "clumpy" , and not azimuthally symmetric about the local vertical at a level which 
can cause many centimeters of error in a baseline measurement. Furthermore, because of incomplete mixing, surface 
measurements are inadequate in estimating this contribution which even at zenith can reach 20 or 30 cm. Ideally, 
this portion of the tropospheric delay should be determined experimentally at each site at the time of the VLBI 
measurements. Often the interferometer data themselves are used to quantify the effect of the water vapor as part of 
the parameter estimation process. 

Independent quantitative estimates of the amount of water vapor along a path through the atmosphere can be made 
by employing "water vapor radiometers" (WVR). These instruments measure the intensity of thermal emission due 
to transitions between rotational energy levels of the water molecule at microwave frequencies (Elgered, 1993). They 
are steerable in both elevation and azimuth, have half-power beam widths of 6°— 9°, and measure sky brightness 
temperatures at several frequencies near 22 and 31 GHz. While there is not a unique correspondence between the sky 
brightness temperature and water vapor content, current data retrieval procedures appear to be capable of accuracies 
as good as several mm of path delay (Keihm and Marsh, 1996). Recent measurements along the lines of sight of VLBI 
observations (Linfield et al, 1996; Teitelbaum et al., 1996) have yielded wet troposphere delays that agree with VLBI 
parameter estimates on the level of a few mm, and give a threefold reduction in residuals. Similar results have been 
obtained by Elgered et al. (1991). Despite optimism in the early years of WVR development, such calibrations have 
not been routinely available for the bulk of VLBI data collected during the past two decades. Because state-of-the-art 
WVR measurements have not been routinely available, VLBI analyses should at the minimum model the neutral 
atmosphere at each station as a two-component effect, with each component being an azimuthally symmetric function 
of the local geodetic elevation angle. 

In an effort to properly account for unmodeled wet troposphere errors, a model of the spatial and temporal spectrum 
of wet troposphere refractivity fluctuations was developed by Treuhaft and Lanyi (1987). This Treuhaft-Lanyi model 
assumes that tropospheric delay errors are dominated by fluctuations in the distribution of water vapor. Furthermore, 
the model assumes that these fluctuations are well described by a Kolmogorov spatial distribution (e.g. Tatarski, 1961) 
that occurs in the bottom 1-2 km of the troposphere and is carried over the site by a constant wind on the order of 
10 m/s. It introduces realistic variations of the troposphere with both time and geometry through this "frozen flow" 
assumption. The model is parametrized in terms of a "rockiness" coefficient, wind direction and speed, and the height 
of the turbulent layer. 

The Treuhaft-Lanyi model is used to generate an a priori observable covariance matrix which is then included in 
the least squares parameter estimation procedure. It provides estimates of correlations of the tropospheric delays 
observed in different parts of the sky at different times. Correlations between phase rate obscrvables are ignored. 
Linfield (1995) demonstrates that they typically do not exceed 10%, and decay to much smaller values on time scales 
of a few minutes. By accounting for spatial variations and not just temporal variations, the Treuhaft-Lanyi technique 
differs from fllters which parametrize the troposphere as a stochastic time- varying parameter (Herring et at, 1990). 
Lack of knowledge of the parameters required by the covariance model presently limits its potential for improving 
VLBI parameter estimates. Some of these parameters, which characterize the strength, extent, and direction of the 
turbulent flow, are starting to be quantified (Naudct, 1995). An important benefit of this technique is that data 
strength is not used to estimate frequent troposphere parameters in cases where the latter are not of primary interest. 
Thus the variances of the estimates of some important parameters {e.g., baselines, source positions) are smaller than 
would otherwise be obtained. This benefit is in addition to the reduction of systematic parameter biases relative to 
estimates which use an inferior observation covariance (data weights). 
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VI. APPLICATIONS AND RESULTS 

The model of VLBI observables which was developed in Sees. Ill— V has been used, with slight variations, by a 

number of research groups to analyze geodetic and astrometric experiments performed with various networks during 
the past two decades {e.g. Kondo et al, 1992; Ma et ai, 1992; Sovers et al., 1993; NEOS, 1994; Johnston et al, 1995). 
The differences among groups generally involve minor variations in tropospheric and tidal modeling. Such analyses 
have presently reached an approximate accuracy level of 1 cm on intercontinental baselines (1 ppb). In this section 
we highlight some of the new insights into a wide variety of geophysics that have resulted from VLBI measurements. 
We intentionally exclude discussion of literature concerning internal structures of extragalactic radio sources and the 
associated astrophysical processes, since it is sufficiently extensive to merit a separate review. Thus, we again focus 
on the accomplishments of purely astrometric and geodetic VLBI. Some of these have already been introduced in the 
sections on model description. We will divide the applications of VLBI into three major categories: reference frames 
(both celestial and terrestrial), structure and dynamics of the Earth, and the orientation of the Earth in the quasi- 
inertial celestial frame of reference. Additional information can be extracted that relates to the Earth's atmosphere, 
relativistic bending, and rudimentary details of source structure. In the near future, data that have accumulated 
over nearly two decades are also expected to permit quantification of the Earth's long-term motion in inertial space 
(galactic rotation). Earlier in this decade, Robertson (1991) presented a survey of VLBI applications that we attempt 
to update here with our brief overview. More detailed reviews, including presentation of recent experiments, analyses, 
and results can be found in Volumes 24 and 25 of Smith and Turcottc (1993), which report the achievements of 
NASA's Crustal Dynamics Project. A general introduction of the applications of many of the new space techniques, 
including radio interferometry, is provided by the book of Lambeck (1988). Much of the astrometric/geodetic data 
that have been accumulated since the 1970s are publicly available from the Crustal Dynamics Data Information 
System at Goddard Space Flight Center. This repository includes both raw data and selected results of analyses, and 
is accessible by computer at http://cddisa.gsfc.nasa.gov/cddis.htinl. 

A. Reference frames 

The accomplishment of VLBI that is of greatest significance is the establishment of reference frames which allow 
quantitative treatment of the dynamics of the Earth and Solar System at unprecedented levels of accuracy. There 
are two such reference systems: celestial and terrestrial. The remainder of this section summarizes determination and 
applications of these two fundamental reference frames, their hundred- fold accuracy improvement since the inception 
of VLBI and related methods, and problems of stability related to these new accuracy levels. 

Astronomical objects have been used for millenia to construct celestial reference systems for measuring the passage 
of time, for navigation, and for investigating Solar System dynamics. To review some of the important milestones in 
the development of celestial reference frames: early astronomers measured the motions of the planets {lit. "wandering 
stars") against the background of "fixed" stars. With improved observational precision, motions of these fixed stars 
became evident. Hipparchus is credited with recognizing precession circa 129 B.C. Proper motions of individual stars 
were observed in 1718 by Hallcy. Nutation of the Earth was discovered by Bradley in 1748. As observing precision 
continued to improve, Herschel and Laplace suggested using extremely distant objects to define astrometric reference 
frames. Such objects reduce or eliminate the effects of the proper motions on reference frame definition. The 1781 
catalog of Messier (Robinson and Muirden, 1979), and the New General Catalog of Dreyer (1888) were important 
steps in identifying these more distant objects. The work of Leavitt (1912) and Hubble (1929) helped to establish the 
extreme distance of what are now classified as extragalactic objects. Present-day Earth-based optical measurements 
have culminated in the latest "fundamental catalog" FK5 (Frickc et at, 1988), and the recent Hipparcos measurements 
from Earth orbit (ESA, 1997) have inaugurated a new era for optical catalogs. 

The historical celestial reference frame, based on measurements at optical wavelengths, has been the unifying 
backdrop for ccnturics-long astronomical observations. It is a coordinate system whose origin is located at the 
barycenter of the Solar System, and which is used to establish locations of objects within the Solar System and 
beyond. Its present accuracy level in the FK5 realization is on the order of 0.1 arcsecond (Fricke et al., 1988). 
Extension of observations into the microwave spectral region by early VLBI measurements already improved this 
accuracy by an order of magnitude. By the late 1990s, yet another order-of-magnitude improvement has brought it 
to 1 mas or better. Unfortunately there are very few objects emitting suflacient flux density to be observable at both 
optical and radio wavelengths. Thus special techniques must be employed in order to connect the VLBI celestial 
reference frame to the historical optical frame (Lestrade et al, 1995). 
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Early determinations of the positions of extragalactic radio sources by Wade (1970) and Cohen and Shaffer (1971) 
yielded angular coordinates accurate at approximately the 1 arcsecond level. Within five years, this accuracy was 
improved by more than an order of magnitude (Wade and Johnston, 1977), and started to rival the best optical 
determinations. Improvement by approximately another two orders of magnitude has been achieved during the past 
two decades for several hundred radio sources. Launch of the astrometric satellite Hipparcos in 1989 (Kovalevsky et 
al., 1995) has likewise extended the accuracy of optical positions of 10^ sources to the milliarcsecond level. Since 
the median optical apparent magnitude of the radio sources is « 18 (well beyond Hipparcos sensitivity), a direct tie 
to the radio reference frame is not possible. A few optical sources (radio stars) also emit weakly at radio frequencies. 
These can be observed in both the optical and radio regions in order to relate the two reference frames (Lestrade et 
al, 1995). 




FIG. 13. Distribution of the 212 bost-obscrvod extragalactic sources comprising the new lAU celestial reference frame (ICRF). 
A conventional equal-area projection is used: declination increases from —90° to +90° (bottom to top), and right ascension 
from —12 to 12 hr (left to right). Note the relatively sparse population at negative declinations. 

During the last two decades, several research groups have used the VLBI technique to catalog the positions of 
extragalactic radio sources (Robertson et al., 1986; Sovers et al., 1988; Ma et al, 1990; Johnston et al., 1995). The 
number of objects whose positions are known at the mas level now exceeds 400. A considerable fraction of these 
continue to be monitored in periodic experiments. Survey campaigns {e.g. Patnaik et al., 1992) with the MERLIN 
and VLA arrays), meanwhile, have accumulated somewhat less precise coordinates of several thousand additional 
radio sources, which form a valuable pool of candidates for inclusion in the higher-precision reference frame. Some 
applications of extragalactic radio reference frames have been to deep space navigation {e.g. Border and Koukos, 1993), 
Earth orientation measurements {e.g. lERS, 1994), geodesy {e.g. Fallon and Dillingcr, 1992), and astrometry {e.g. 
Bartel, 1985; Treuhaft and Lowe, 1991; Lebach et al., 1995). In addition to the intrinsic scientific interest in the stability 
of dynamical systems, these diverse applications require accurate and stable positions of the objects composing the 
reference frames. The International Astronomical Union has now adopted a new fundamental celestial reference frame 
(the International Celestial Reference Frame, ICRF) that is based on the angular coordinates of 212 radio sources 
(lAU, 1996; Ma et al., 1997). Figure 13 illustrates the distribution of the defining sources over the sky. This is the 
first time that the fundamental celestial coordinate system is no longer based on observations at visible wavelengths. 
The ICRF is approximately 100 times more accurate than the FK5 catalog, the present realization of the fundamental 
optical celestial frame. 

Limitations of reference frames that are based on extragalactic radio sources have been recognized for some time. 
While estimates of the proper motion of quasars presently give null results at the approximate level of 50 prad/yr 
(Eubanks et al., 1996), the extended structure of the emitting objects at the nrad level is ubiquitous. Structure that is 
constant in time can be handled by careful definition of a fiducial point for each source. Unfortunately, many sources 
show considerable time variability, and thus require constant monitoring to ensure stability of the reference frame. 
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Such variability can be considered to belong in the category of systematic errors. Specialized imaging experiments 
are being carried out in order to correct for source structure and its variation (Fey et al., 1996). Imaging is also 
feasible, however, via reanalysis of older data collected primarily for geodetic purposes (e.g. Chariot, 1994; Finer and 
Kingham, 1997). Figure 14 shows such results for the source 3C 273 in 1986. It can be seen that the central core 
and neighboring components ("jet") extend over several mas, making it impossible to establish a fiducial mark at 
the sub-milliarcsecond level without detailed study. Similar images at other epochs show that the components of 
the jet move rapidly, and new components are ejected sporadically. After this structure and its time variation was 
characterized in the 1980s, it was dropped from most observing schedules. The vast majority of the defining sources 
of the ICRF, however, are much more point-like. 



3C 273 




Right ascension 



FIG. 14. Typical maps of the S- and X-band intensity distribution of the source 3C 273, generated from geodetic VLBI data 
for the epoch 1986.27. 

Another source of instability of the extragalactic reference frame originates in the motion of the Solar System due 

to the rotation of our Galaxy relative to the extragalactic radio sources. Detection of the acceleration of the SSB 
is now on the borderline of possibility. As mentioned in Sec. III.E, the motion of the Solar System in inertial space 
makes its presence felt if observational data extend over a sufficiently long time span. To the extent that the Earth 
moves along with the uniform rotation of our Galaxy and the distance to the Galactic Center is known, this would 
amount to an indirect measurement of the period of galactic rotation. Such results could also make a contribution 
to our understanding of cosmology from comparisons to motion of the SSB with respect to the cosmic background 
radiation. A related topic is the proposed use of VLBI measTirements to determine the distance to the Galactic center 
(Reid, 1993) which would be of importance in establishing a more reliable cosmic distance scale. 

The terrestrial reference frame has its origin at the Earth's center of mass and establishes fiducial points on the 
surface of the Earth. Prior to the advent of space geodetic experiments, there was no global connection of national 
geodetic grids at a level better than ftil meter. Doppler tracking of artificial Earth satellites improved this situation, 
and the advent of laser ranging, VLBI, and global positioning satellite tracking has further increased the accuracy 
of global reference frames to approach 1 cm in the late 1990s. In the arena of terrestrial reference frames, VLBI 
measurements have contributed significantly to the establishment of more than 100 fiducial points fairly uniformly 
distributed over the Earth's surface, which collectively comprise the ITRF: International Terrestrial Reference Frame 
(Boucher et al., 1996). Figure 15 shows the global distribution of VLBI stations that contribute to the ITRF, while 
Fig. 19 of Sec. VLB gives an example of such results. The VLBI technique has not been the only contributor to an 
improved terrestrial reference frame. Laser ranging to satellites and the recent proliferation of GPS measurements 
have begun to play a dominant role in forming a unified terrestrial coordinate system that is accurate at the centimeter 
level. This provides a global reference for the national geodetic networks which have been traditionally established 
by conventional triangulation methods; these networks are being densified by GPS measurements. The result is a 
globally coherent reference system which provides a backdrop for numerous studies of the Earth's behavior. Some 
limitations to its stability are discussed in the next section. 
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FIG. 15. Distribution of 104 VLBI observing stations that contribute to the International Terrestrial Reference Frame (ITRF). 
An equal-area projection of the Earth's surface is used. Note the heavy concentration in the northern hemisphere. 

It must be noted that both current celestial and terrestrial reference frames exhibit a pronounced northern bias. 
Due to the dominance of observing stations in the northern hemisphere, there has been nonuniform coverage of both 
the Earth's surface and the celestial sphere. For example, fewer than one fifth of the 50 best-determined ITRF sites 
are in the southern hemisphere, and fewer than one third of the 212 sources that are included in the new ICRF are 
at negative declinations (Figs. 13 and 15). It is expected that future development of observatories in the southern 
hemisphere will remedy this situation to the extent possible with the limited land area available south of the Earth's 
equator. 

B. Earth orientation and structure 

The motion of the Earth as a whole, as well as its surface and bulk composition and structure, are also areas in which 
VLBI measurements have made significant contributions. There are basically two classes of Earth-related applications. 
They pertain either to orientation of the structure as a whole, or to motions of its crust. Both categories have varied 
root causes: Solar System dynamics and internal, oceanic, and atmospheric processes on the Earth. Manifestations 
of most of these complicated motions are for the most part far removed from ordinary human experience. One 
exception may be the irregular slowing of the Earth's rotation, which necessitates the introduction of occasional "leap 
seconds" in order to keep synchrony between the Earth's rotation and atomic time. Many of the remaining effects 
have amplitudes that are only on the order of a meter or less at the Earth's surface, but their characterization is of 
crucial importance in building a coherent picture of our planet and its environment. 

Precise specification of the orientation of the Earth is of practical importance, for example, in navigating spacecraft 
to other planets. The seemingly minuscule errors of parts per billion are magnified into errors of kilometers near 
Jupiter. Predictions of future Earth rotation values (particularly UTl) lose accuracy very rapidly for times past a 
few weeks after the date of measurement and extrapolation. Continuous monitoring of Earth rotation is therefore 
imperative for applications that require high-quality real-time UTPM. Substantial help in this regard comes from a 
perhaps unexpected source: the global distribution of atmospheric winds or more precisely total atmospheric angular 
momentum (AAM), which is widely monitored in real time for weather forecasting. 

Measurements of the rotation of the Earth and the wandering of its spin axis relative to the crust contain a great 
deal of information related to diverse physical processes. While to a good approximation the rotation rate is constant, 
it has slowed considerably over geological time scales. The length of a day was only 18 hours 900 million years 
ago as determined from recent analyses of tidal sediments (Sonett et al., 1996). A number of phenomena come into 
play in determining the rotation rate, ranging from the lengthening of the Earth-Moon distance, to frictional and 
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electromagnetic forces between the Earth's core and mantle, to frictional forces between winds and ocean water and 
the surface. Presently the rate of deceleration is much larger than the billion-year average, on the order of 2 ms 
per day, but it is known to have reversed sign within the last hundred years (Archinal, 1992b). The spin axis has a 
quasi-periodic motion that describes an approximate circle of w 20 m diameter (see Fig. 16) in w 300 days whose 
center also wanders by tens of meters over decades; the size of the circle fluctuates with a period of approximately 7 
years. In addition to motions with respect to the Earth's crust, the spin axis also exhibits motion in inertial space 
(nutation). The two nutation angles arc usually included along with the three quantities describing Earth rotation 
to form a set of five "Earth orientation" parameters. The radio interferometric technique, along with other space 
geodetic experiments, has brought about orders-of-magnitude improvements in measurements of small irregularities 
in these motions, thereby opening prospects for detecting numerous periodic and aperiodic processes. We summarize 
the present accuracy of Earth orientation determinations, and present examples of contributions of the results to two 
diverse fields of geophysics. 

With a network of several stations that are reasonably uniformly distributed over the surface of the Earth, and 
separated by baselines on the order of an Earth radius, it is possible to infer the Earth orientation parameters to an 
accuracy level that is better than 1 mas (a few nanoradians, or parts per billion). To illustrate VLBI measurements 
of Earth orientation. Figs. 16 and 17 show the variation of the position of the spin axis relative to the crust, and the 
rate of rotation, during a part of the 1980s and 1990s. As seen in Fig. 16 the spin axis follows a nearly circular path, 
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FIG. 16. Observed crust-fixed position of the Earth's spin axis during the years 1986-1993. The origin is at the conventional 
pole of 1903, and the upward arrow indicates the approximate present-day average pole position. 

with a diameter that oscillates between a few meters and 20 m with a period of approximately seven years. The 
deceleration of the Earth's rotation rate is seen in the bottom part of Fig. 17 as a lag behind atomic time (TAI). 
The resulting excess length of day is seen in the top part of the figure: it exhibits a rich spectrum which is just now 
beginning to be understood in detail. The dominant oscillations are caused by the annual interchange of atmospheric 
angular momentum with the solid Earth. 

Experimental determinations of the Earth's rotation rate show that it is not constant at the level of a few parts per 
billion. There are annual and semiannual cycles in the length of day, both with amplitudes of approximately 0.3 ms 
and minima in January and July (Eubanks, 1993). Since the early 1980s it has been known that most of this variability 
is well correlated with the total AAM (Morgan et al., 1985). Refinement of this correlation was made possible by the 
superior precision of VLBI relative to classical Earth orientation measurements using optical techniques. Variations 
in the atmospheric angular momentum are transferred (via friction and forces on mountain ranges) to variations of 
the angular momentum of the solid Earth, conserving total angular momentum. This discovery has made substantial 
improvements in our ability to perform short-range forecasts of Earth rotation. Global weather measurements play 
a significant role in this process, and their rapid availability is crucial in improving forecast quality (Freedman et 
al., 1994). As with most global geophysical measurements, the spatial and temporal density is inadequate at some 
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level. Relatively few weather stations operate in the southern hemisphere, and the strong winds at high altitudes 
(pressures of 100 to 1 mbar) are significant (Rosen and Salstein, 1985). In addition to the dominant 0.6-ms effect of 
atmospheric angular momentum, a number of smaller meteorological and tidal effects contribute to irregularities in 
Earth rotation. Before considering some selected examples of these, we discuss the remaining components of Earth 
orientation; the two nutation angles that specify Earth orientation in an inertial frame. 
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FIG. 17. Observed rotation rate of the Earth during the years 1980-1996. The bottom plot shows the discrepancy of time 
measured by atomic clocks and by the Earth's rotation, while the top plot shows the increase of a few ms per day in the length 
of the day. The annual variations are caused by wind-surface friction. 



In response to torques exerted by the Sun and Moon, the Earth's spin axis oscillates in fundamental modes dis- 
tributed over a wide range of frequencies. The amplitudes of these motions can be calculated fairly precisely based on 
classical dynamics. At the discrimination level of modern measurement techniques, however, the theoretical models 
are not adequate. The relatively poorly understood internal structure of the Earth plays a large role in some of 
these motions, and empirical measurements have jiist begun to shed light on the underlying mechanisms. In the 
1970s it became increasingly clear that the then-standard Woolard nutation model was insufficiently accurate for 
analyses of modern data. Based on the theoretical work of Gilbert and Dziewonski (1975) and Kinoshita (1977), 
Wahr developed a new nutation scries (1979; 1981), which was adopted by the International Astronomical Union in 
1980 (Scidclmann, 1982). These two series describe the time variation of the nutations in longitude and obliquity by 
superposing oscillations at a total of 106 frequencies, with amplitudes of 0.1 mas precision. Only a few years after the 
adoption of the 1980 lAU model, VLBI measurements revealed significant errors, predominantly at semiannual, an- 
nual, and 18.6-yr periods. Herring et al. (1986) initially pointed out these discrepancies, and gave improved estimates 
of a number of amplitudes from VLBI analyses. Gwinn et al. (1986) subsequently interpreted them as manifestations 
of irregularities in the figure of the Earth's inner core. 

As an example of the ability and evolution of VLBI techniques in the measurement of the Earth's spin axis in 
inertial space, Fig. 18 shows plots of the two nutation angles Ail) and Ae during the history of regular astrometric 
VLBI measurements. The plots show corrections to the standard (1980 lAU) theoretical model of nutation. It is 
apparent that there are highly regular deficiencies at the 10-mas level. Yearly and semiannual frequencies are clearly 
visible, along with long-term (linear and 18.6-y period) discrepancies of considerable size. 

The dominant nutation in longitude, with a period of 18.6 years, is difficult to separate from precession (26,000- 
yr period) with short data spans. Only recently has the history of VLBI observations reached a time span of one 
18.6-yr period, permitting this decoupling. Determinations of the nutation amplitudes and precession at the levels of 
a milliarcsccond and a few tenths of mas/yr, respectively, have been made by Herring (1988), Chariot et al. (1995), 
Walter and Ma (1994), and Walter and Sovers (1996). The value of the precession constant which was determined 
from optical measurements in the 1970s (Pricke, 1977) has been found to be in error by 3 mas/yr. Most recently, it has 
been found that torques from bodies other than the Sun and Moon can contribute to nutation of the Earth at a level 
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of 0.1 mas (Kinoshita and Souchay, 1990). These have not yet been directly detected, but evidence is accumulating 
that accounting for their cumulative effect improves the fit of the VLBI model to experiment. Another revision of the 
official lAU nutation series is expected to include several hundred additional terms, and to correct many amplitudes 
on the basis of empirical VLBI determinations. 
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FIG. 18. Discrepancies of measured nutation angles and the 1980 lAU models during the years 1984-1996. Both the obhquity 
Ae and longitude A^sine are shown. Short-term (annual and semiannual) and long-term (linear and 18. 6- year) discrepancies 
are evident. 

Gravitational forces from the Solar System, predominantly the Sun and Moon, cause a rich spectrum in the tidal 
response of the solid Earth and its oceans and atmosphere. Direct motions of the Earth's crust (with amplitudes 
on the order of a meter) arc quite obvious in VLBI analyses, but smaller secondary effects on station positions and 
Earth orientation also cause measurable effects, and need to be considered. One of the first estimates of the solid 
Earth's tidal response was the determination of Love numbers from VLBI by Herring et al. (1983): h = 0.62 ± 0.01, 
I = 0.11 ± 0.03, and a phase lag = 1° ± 1°. More recent work (Haas and Schuh. 1996; Gipson, 1996) shows that the 
response factors at numerous tidal frequencies can be reliably obtained from analyses of VLBI measurements. 

Perhaps the most important of the secondary tidal effects is the modification of Earth orientation by mass redistri- 
bution in the oceans by global tides. This affects both the rotation rate (UTl) and orientation of the spin axis (PM). 
Brosche et al. (1989) made the initial attempts to model these "fast UTPM" variation amplitudes. In the early 1990s 
several independent determinations by VLBI and SLR of the Earth orientation amplitudes induced by 8 ocean tidal 
components gave results in good agreement with each other (Sovcrs et al., 1993; Herring and Dong, 1994; Watkins et 
al, 1994). Initially their agreement with purely theoretical values based on global tide models was not good, but it 
is improving (Gross, 1993; Ray et al, 1994; Chao et al, 1996; Gipson, 1996), as the global tide model is refined by 
TOPEX/Poscidon results. Agreement is now at the level of 10/xas. 

In addition to affecting the orientation of the whole Earth, global redistribution of mass in the oceans by tidal action 
also causes local motions of stations on the Earth's crust. These "ocean loading" amplitudes have been incorporated 
in VLBI models since the 1980s (Scherneck, 1983), and can amount to motions of several cm. Although their small 
size places them close to the current resolution limit of VLBI, long data spans are capable of yielding significant direct 
determination of their amplitudes (Sovers, 1994). By analogy with oscillations in the distribution of ocean water, 
variations in the atmospheric pressure at an observing site can also cause long-term vertical motions with amplitudes 
of many mm. This "atmospheric loading" has been detected in VLBI analyses (Rabbel and Schuh, 1986; Manabe et 
al, 1991; van Dam et al, 1994). 

Another minor systematic effect which has only very recently been identified is post-glacial reboimd: the slow 
relaxation of the Earth's crust in response to melting of the thick ice sheets which covered it during the last glacia- 
tion. The dominant motion is vertical, and theoretical models produce estimates of several mm/year (Peltier, 1995). 
Unfortunately, determination of the vertical component of station motion is notoriously weak relative to the two 
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horizontal components, because of the strong correlation of the vertical with the zenith tropospheric delay. Neverthe- 
less, recent empirical estimates for a number of VLBl observing stations with sufficiently long observation histories 
(Argus, 1996; Ryan et al., 1997) have produced results that are in fair agreement with theoretical models (Mitrovica 
et al, 1993). 

A limitation to the stability of the terrestrial reference frame, which is also of great interest for its own sake, is 
the motion of the tectonic plates that constitute the Earth's solid surface. Since this possibility was proposed by 
Wegener in the early part of the 20th century, the existence of tectonic motion has been confirmed for time scales 
of millions of years. This was achieved by detailed cataloging of magnetization properties of crust near active plate 
boundaries (Minster and Jordan, 1978; DeMets et al., 1990). As the time base for space geodetic station motions 
became sufficiently long in the 1980s, the astonishing result emerged that present-day (decade scale) plate motions are 
very similar to those inferred from paleomagnetic data (Myr scale). This attests to the relatively smooth character 
of the forces driving plate tectonics. At present, the rates of motion of more than a dozen tectonic plates that have 
been occupied by VLB! stations are known to within fractions of mm/year (Fallon and Dillinger, 1992; Argus and 
Gordon, 1996; Boucher et al, 1996). 
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FIG. 19. Measured baseline between Alaska and Hawaii during 514 VLBI observing sessions, 1984-1997. L is the baseline 
length in mm, minus a constant offset of 4728114 meters. The weighted RMS scatter about the linear trend of —46 mm/yr is 
7 mm. 

As an example of the manifestation of tectonic station motion, Fig. 19 shows results of VLBI measurements of the 
length of the baseline connecting stations at Gilmorc Creek, Alaska and Kokce Park, Hawaii (island of Kauai) during 
the years 1984-1997. This figure is derived from results of Goddard Space Flight Center 1997 data analyses, archived 
at the World Wide Web site http://lupus.gsfc.nasa.gov/vlbi.htinl. Several points of interest are illustrated by 
the time history of this particular baseline. First, the North American (Alaska) and Pacific (Hawaii) plates arc one of 
the fastest moving pairs of tectonic plates, converging at nearly 5 cm/yr. Second, there is nearly perfect agreement 
of the measured rate of change (—46 mm/yr) with the Nuvel model rate of —45 mm/yr. Third, the weighted root- 
mcan-square scatter about a perfectly linear trend is approximately 7 mm, which is a typical result for a baseline 
of this length («4700 km), and amounts to 1.5 ppb. Fourth, the evolution of the VLBI technique is illustrated by 
the increased frequency of experiments and the decrease in baseline length uncertainties during the 1980s and 1990s. 
Fifth, even the current theoretical models have some systematic errors, which may be manifested in some of the 
departures from linearity, e.g., in 1994-96. Finally, the general stability of the experimental techniques is shown by 
the smooth transitions at the epochs of both a large natural and a man-made perturbation at Kauai: the hurricane 
damage in late 1992 and the transition to a new antenna, 40 meters distant, during the summer of 1993 are hardly 
noticeable in Fig. 19. 

Proximity of VLBI stations to the epicenters of strong earthquakes, and to active faults (especially the San Andreas 
fault in California) is starting to contribute to improved understanding of the forces that are involved in these motions. 
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Several pairs of measurements of station coordinates before and after substantial earthquakes now exist (Clark et 
al., 1990). They confirm permanent, essentially instantaneous displacements of several cm, but the resolution is not 
sufficient to detect any relaxation mechanisms in the subsequent hours and days. The relatively dense network of 
stations on both sides of the San Andreas fault show that the magnitude of horizontal motion associated with the fault 
movement decreases with increasing distance from the fault (Ward, 1988). Two reviews (DeMets, 1995; Larson, 1995) 
contain more detailed discussion of recent results related to global tectonics and crustal deformation. A recently 
developed technique of active interferometry. synthetic aperture radar (SAR), shows promise of accurate and timely 
determination of regional deformations associated with tectonic and earthquake activity (Zebker et al., 1994). 

C. Troposphere and ionosphere 

Atmospheric effects have long been the bane of Earth-based astronomy. They continue to play that role for 
interferometry at radio wavelengths. In contrast to many of the more deterministic parts of the VLBI model, the 
effects of the atmosphere are not easily predicted or externally quantified. The experiments themselves must therefore 
serve to characterize atmospheric behavior. In addition to the global effects of the atmosphere (AAM) that were 
already mentioned in connection with the above discussion of Earth rotation, we consider two additional aspects 
here. At best, empirical estimation can yield time series of accurate atmospheric parameters, whose applicability is 
unfortunately limited to the times and locations of the experiments. 

Determination of ionospheric electron content is possible whenever the VLBI measurements are performed at two 
widely separated frequencies. The details of this dual-frequency calibration process can be used to determine the 
integrated electron density along the incoming ray path. Limited sampling in both time and space, however, limits 
the utility of VLBI as a method for studying the global ionosphere. Satellite techniques give a much denser data 
set, and routine generation of ionospheric maps is currently being implemented using GPS measurements both with 
ground-based receivers and other satellites in low Earth orbit (Wilson et al., 1995; Ho et al., 1996). 

Estimating the zenith tropospheric delay at frequent time intervals for each VLBI station during data analysis 
permits monitoring the time variation of water vapor content above that station. These intervals can even be reduced 
to essentially the intervals between observations if a stochastic filtering scheme is used {e.g. Herring et al, 1990). 
If the observing schedule has ensured a good spatial distribution of measurements (including low elevation angles), 
then the uncertainties of the derived zenith wet delays are typically on the order of a few mm. Comparisons of water 
vapor radiometer and VLBI results have demonstrated the ability to track the atmospheric water vapor content quite 
faithfully (Elgered et al., 1991; Teitelbaum et al., 1996). For the highest accuracy, it is important to account for 
both the temporal and spatial correlations of fluctuations of the tropospheric delay. Stochastic estimation schemes 
can account for temporal correlations. To describe atmospheric spatial asymmetry, horizontal gradient parameters 
may be estimated. An analysis scheme utilizing a theoretical model of atmospheric turbulence which also contains 
parameters to quantify the "dumpiness" of the atmosphere (Treuhaft and Lanyi, 1987; Treuhaft and Lowe, 1991) 
accounts for correlations between VLBI delays as a function of both time and space. It also gives a better description 
of the influence of both the water vapor content and wind distribution surrounding a site (Naudet, 1995). 

D. Relativity 

At the precision level of intercontinental VLBI measurements, the incoming signals are affected by their passage 
through the varying gravitational potential within the Solar System. Since the positions and velocities (ephemerides) 
of all the major gravitating bodies are known to a high degree of accuracy, analyses of the experiments provide a 
means of testing alternative relativistic theories. Shapiro (1964) first suggested that timing experiments on radio waves 
propagating through the Solar System can be used to measure the parameter 7ppjj of parametrized post-Newtonian 
relativity theory. These original measurements of radar reflections from the inner planets (Shapiro, 1967) gave the 
result 7ppN — 0.9 ± 0.2. Determinations of 7ppN from geodetic VLBI data spanning more than a decade gives 7ppN 
= 1.000±0.002 (Robertson and Carter, 1984; Robertson et al., 1991). This uncertainty, = 0.002, is competitive 
with the results of specifically designed experiments of much shorter duration. An example of the latter is a VLBI 
experiment which also observes at a third frequency (K band, >12 GHz) in addition to the standard S-band and 
X-band wavelengths (Lebach et al., 1995). This yields 7ppj^ = 0.9996±0.0017, a somewhat improved despite 
an order of magnitude fewer data. Experiments scheduling numerous observations in the vicinity of the Sun can 
achieve equivalent results with just a fraction of the data, because of the extremely steep variation of gravitational 
delay with closest approach to the Sun. The higher K-band frequency is less sensitive to the corrupting influence of 
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plasma-induced phase scintillation. The estimated value of unity for 7ppN (within one standard deviation in all cases) 
indicates agreement with Einstein's general relativity. 

Unfortunately none of the other parameters of post-Newtonian relativity theory are accessible to measurement 
with VLBI experiments. Nonlinearity in the superposition of gravity [3 and the time dependence of the universal 
gravitational constant G/G cannot be determined from VLBI, but limits are being placed on their magnitudes by 
analysis of lunar laser ranging data (Williams et al., 1996). Initial searches for VLBI evidence of the passage of 
gravitational waves have placed loose boimds on their existence (Pyne et ai, 1996; Gwinn et ai, 1997). We also note 
that VLBI is weakly sensitive to the difference between Newtonian and Lorentzian aberration [oc l/4(i;/c)^], which is 
a 1-cm effect for long baselines and is implicitly included in the VLBI model. 

VII. PROGNOSIS FOR FUTURE MODEL IMPROVEMENTS 

Data analysis has provided continuous feedback to experiments by characterizing previously unknown or poorly 
quantified aspects of the model (e.g., nutation, tidal variations, troposphere). Such interplay is certain to continue 
in the future, as both experimental and theoretical techniques are refined to eliminate remaining systematic errors. 
This section gives condensed summaries of areas which will require close attention in the future if improvements of 
the current VLBI model are to permit achieving true part-per-billion accuracy. 

Both special and general relativistic aspects of the current model could be improved. Second-order general rel- 
ativistic effects have not yet been thoroughly investigated, but probably do not contribute at the picosecond level. 
Theoretical studies such as that of Turyshev (1996) may make a valuable contribution here. If variations of the 
gravitational potential along the path of the baseline through the Earth are taken into account in calculating the 
proper distance, this correction was estimated by Thomas (1991) to amount to 2 mm for a 10,000 km baseline. Sim- 
ilarly, variations in the gravitational potential at the; station clocks are only approximately accounted for by means 
of Eq. (3.31). Concerning special relativistic aspects, the Fairhead and Bretagnon (1990) extension of the work of 
Moyer (1981) on the "time ephemeris" produces higher-order terms that contribute to TDB— TDT at the ^s level. 

Galactic effects may soon emerge above the detection threshold. The rotational motion of the Galaxy produces 
aberrational effects which change by 20 prad/yr. This will need to be taken into account for observations spanning 
more than two decades. 

The rich tidal spectrum of the solid Earth and its oceans will be a source of model refinements for some time to come. 
Direct contributions of the planets to solid Earth tidal displacements can reach the millimeter level. In addition to the 
eight frequencies considered in the model described in Sec. III.D.l.a, short-period variations of UTPM have additional 
components (Seller and Wiinsch, 1995: Gipson, 1996). Those which are significant at the mm level will emerge as 
data analyses are refined. Empirical estimates of ocean loading amplitudes for several IRIS stations (Sovers, 1994) 
indicate that the best theoretically derived amplitudes may be in error by several mm. Future refinements in data 
analyses, and improved global ocean models from the recent TOPEX/Poseidon mission (Le Provost et at, 1995; Fu 
et al., 1994) may improve the accuracy of the theoretical ocean loading model to the mm level. Resonance with the 
Earth's free core nutation may modify some of the amplitude corrections at nearly diurnal frequencies by « 1 mm. 
Ocean tides cause motion of the center of mass of the solid Earth due to motion of the center of mass of the oceans 
(Brosche and Wiinsch, 1993). The amplitude of this displacement can be as large as 1 cm at the usual diurnal and 
semidiurnal tidal frequencies. Non-conservation of the total angular momentum of the Earth via interchange with 
the Sun and Moon is also not included in the tidal models (Brosche and Seller, 1996). The effects of both of these 
refinements on VLBI observations must be assessed. The retarded tidal potential effect mentioned in Sec. III. C. 2. a 
can be as large as several tenths of mm. Thus, for correct modeling at the mm level, the light travel time should be 
accounted for. 

Non-point like flux distributions of the radio sources are being efficiently mapped by experimenters using the 
VLB A (Fey et al., 1996). Their results are expected to provide improved models for astrometric and geodetic VLBI. 
Meanwhile, estimates of parameters for simplified structural models may improve data analyses. 

There are short-period deficiencies in the present International Astronomical Union models for the orientation of 
the Earth in space that may be as large as 1 to 2 milliarcseconds, and longer-term deficiencies on the order of 1 mas 
per year (3 cm at one Earth radius) . VLBI measurements made during the past decade indicate the need for revisions 
of this order of the annual nutation terms and the precession constant (Eubanks et al., 1985; Herring et al., 1986). 
The 18.6-year components of the lAU nutation series are also in error, and present data spans are just approaching 
durations long enough to separate them from precession. Options to improve the nutation model were discussed 
in Sec. III.D.2. Any of these constitute a provisionally improved model, especially for the annual and semiannual 
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nutations, until the lAU scries is officially revised. Future refinements of the equation of the equinoxes (Eq. 3.86) will 
probably lead to changes on the order of tens of i^s in the hour angle. 

The large dish antennas that are used to collect signals from extragalactic sources are susceptible to various insta- 
bilities at the level of many mm. Gravity loading may cause systematic variations in the position of the reference point 
of a large antenna that are as large as 1 cm in the local vertical direction. Such systematic errors and their dependence 
on antenna orientation and temperature may be modeled (Clark and Thomsen, 1988; Jacobs and Rius, 1989). The 
geometric structure of each antenna, as well as its alignment with respect to local site features, should be carefully 
checked against design specifications. For example, hour angle misalignment on the order of 1 arc minute can cause 
1 mm delay effects for HA-Dec antennas with 7-m axis offsets. Thermal expansion of the portion of an antenna 
above the reference point may induce delay signatures that are several mm pcak-to-peak for a typical VLBI antenna 
(Nothnagel et al., 1995). It is thus imperative to model this effect for achievement of the highest accuracy. 

The limits of validity of the dual-frequency calibration procedure for ionospheric effects need to be carefully estab- 
lished, in conjunction with consideration of plasma effects for ray paths near the Sun. In addition, corrections for the 
gyrofrequency effect may reach a millimeter. 

New techniques for characterizing the atmosphere are expected to allow more realistic modeling of the tropospheric 
delay than the simple spherical-shell model underlying all the results of Sec. V.B. When comprehensive atmospheric 
data from a region surrounding each observing site are available, present computer speeds should permit estimating 
the tropospheric delay by means of a complete ray-tracing solution for every observation. Meanwhile, improvements 
in tropospheric mapping can be sought by modeling variations of the temperature vs. altitude profile as a function 
of season, latitude, altitude, and diurnal cycle. Efforts are also under way [e.g. MacMillan and Ma, 1997; Chen and 
Herring, 1997), to model azimuthal gradients in the troposphere. Persistent equator-to-pole gradients in pressure, 
temperature, humidity, and tropopause height suggest that a priori modeling of North-South gradients may be 
beneficial. East- West gradients, which are probably dominated by weather systems passing over a site, are likely to 
be more difficult to model without extensive weather data. 

By the late 1990s astrometric and geodetic radio interferometry has largely fulfilled the promise seen for the tech- 
nique at its inception in the late 1960s. VLBI has yielded a new nearly inertial celestial reference frame that is 
accurate at a level of nanoradians, achieved point positioning on the Earth at the centimeter level, and produced 
the capability of determining the instantaneous orientation of Earth in space at similar accuracy levels. The VLBI 
technique has intersected many fields of physics, and its evolution has involved both routine and unexpected com- 
ponents. In addition to improved quantification of known physics, there was need to take into account a number 
of unexpected contributions, which led to enhanced understanding of the behavior of the Earth. The field is still 
developing vigorously both on the experimental and theoretical fronts. It can be hoped that progress will not be 
unduly impeded by numerous processes effective at the millimeter level. The future holds the promise of exciting and 
unexpected results characterizing the behavior of the Earth, the structure and emission mechanisms of extragalactic 
radio sources, and the motion the Earth in the universe. 
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TABLE I. Maximum Magnitudes and Present Uncertainties of Portions of the VLBI Delay Model (mm). 





IxyTp-viTTmiTi Mf^lmr 
ivxcLA.1111 Lliii ucicty 












11 n ppTt ?i i n tv 

LtlX^^CX \jCXiLLL\jy 


RA^PT TNP HFOA/TFTRV 










fi V 1 






hjcll ViLL ui uitdi iiiutujii 


U A 1\J 




1 


Vjri dVlLdLKJllcll Llcidy 


9 V 

^ A ±U 




9 


STATION POSITIONS 








-LcL-tUillL- X11(JI/1U11 


1 nn 




1 

i. 


TidSil motion 


ouu 







iNOli liCldi IIIOLIUU 


ou 




c; 











UTPM 


2 X 10** 




2 


Nutation /precession 


3 X 10^ 




3 


SOURCE STRUCTURE 


50 




10 


ANTENNA STRUCTURE 


10* 




10 


INSTRUMENTATION 


3 X 10^ 




5 


ATMOSPHERE 








Ionosphere 


10^ 




1 


Troposphere 


2 X 10* 




20 


TABLE II. Tectonic Plate Rotation Velocities: NNR-Nuvel-IA Model 
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TABLE III. Frequency Dependent Sohd Earth Tide Parameters. 
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TABLE IV. Lunar Node Companions to Ocean Tides. 
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TABLE W Ocean Tidally Iiuluccd Pt-riodic Xariatious in I'olar Motion (Sov(-rs <i al.. 19!):5). 
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TABLE \T. Solar System Wlocity Relative to \ arious Standards of R(>st.. 



Reference 


Velocity (km/s) 


I (deg) 


b (deg) 


Local standard of rest*^ 


20 ± 1 


57 


23 


Galactic center''''' 


240 ± 5 


90 





Local Group'^ 


308 ± 23 


105 


-7 


Cosmic microwave background'^ 


370 ± 3 


264 
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Kerr and Lynden-Bell (1986) 
Fich et al. (1989) 
" Yahil et al. (1977) 
Kogut et al. (1993) 
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TABLE VII. Plasma Effects. 



Plasma p{e/rrC') 


fp (kHz) 


K/t's) 


(J^f/J^x) 


Earth 10'^ 
Interplanetary 10^° 
Interstellar lO'^ 


8980 
898 
28 


4 X 10"^ 

4 X 10"" 
1 X 10-5 


10"^ 

10-* 
3 X 10"® 


TABLE VIII. Electron Gyrofrequency Effects. 


Magnetic field B (gauss) 


(kHz) 






Earth 0.3 
Interplanetary 10~^ 
Interstellar 10"® 
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0.028 
0.003 


4 X lO"" 
1.2 X 10"® 
1.2 X 10"^ 


10-" 
3 X 10-'^ 
3 X 10"^° 


TABLE IX. Sensitivity of Tropospheric Delay to Lanyi Mapping Function Parameters. 


Parameter 


Standard value 




Sensitivity (6°) 


To 
W 
hi 
h2 


292 K 
6.8165 K/km 
1.25 km 
12.2 km 




-7 mm/K 
20 mm/K/km 
—20 mm/km 
5 mm/km 
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